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PREFACE 


A Kalman filter was designed to yield optimal estimates 
of geophysical parameters from (observed) VLBI group delay 
data. The geophysical parameters are the polar motion 
components, adjustments to nutation in obliquity and longi- 
tude, and a change in the length of day parameter. VLBI 
clock parameters and atmospheric zenith delay parameters 
were simultaneously estimated with the geophysical param- 
eters. The right ascension and declination of the VLBI 
(radio) sources, and VLBI site positions were kept at a 
priori values. The primary source of VLBI group delay 
observations was the 1984-85 IRIS (International Radio 
Interferometric Surveying) data set consisting of day long 
VLBI observing sessions (at five day intervals) . The Kalman 
filter limits subjectivity in VLBI data analysis. 

The scientific objective of Kalman filtering the VLBI 
group delay data was to examine the relationship between 
earthquakes and associated mass movements, and the Chandler 
wobble of the earth. The Kalman filter produces precise 
estimates of polar motion over time intervals of one day or 
less. The estimated polar motion results were deconvolved 
to produce excitation functions versus time. Changes in the 
earth's inertia tensor due to earthquakes and associated 
mass movements may be seen in the observed excitation 
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functions. The excitation signatures of the preceding 
mass movements may appear as Heaviside step functions or 
other manifestations. Such excitation studies have been 
made in the past using less precise polar motion data; one 
objective of this research was to characterize excitation 
signatures associated with earthquakes more exactly using 
contemporary data. 

Several great earthquakes occurred during 1984-86, all 
with magnitudes (M s ) approaching 8.0 or above. Excita- 
tion functions were calculated at five day intervals around 
the time of the 1985 great Chilean and 1985 great Mexican 
earthquakes. The excitation results were produced for a 
45 day window surrounding the epoch of each main event. 

The observed excitations were compared with theoretical 
excitation (elastic dislocation) values calculated by Chao 
and Gross (1987). 

Changes (ramps) in excitation are strongly correlated 
with the time of the 1985 great Mexican earthquake, but 
the changes in excitation are not statistically significant. 
The excitation formal errors are typically 0.5 milliarc- 
seconds. The observed excitation magnitude and direction 
are very close to being in accord with the theoretical 
estimates. 

Two steps are observed in the x-component of excitation 
twenty days prior to the 1985 great Chilean earthquake and 
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roughly twenty days after. The step magnitudes are 1.5 mas 
and 0.8 mas, respectively. The net excitation is 2.4 mas 
in a direction 7.1 degrees East longitude. This is signi- 
ficant. The observed excitation is thirteen times larger 
than the theoretical value. The theoretical excitation 
direction is 110 degrees East longitude. These Chilean 
results are similar to the findings of Chao and Gross 
(1985) for the great 1977 Sumba earthquake. 

Finally, IRIS polar motion data are deconvolved near the 
time of the great 1986 Taiwan earthquake. A step of 3.8 
+/- 0.5 mas in the (observed) x-component of excitation 
begins directly after the Taiwan main shock. Such a step 
is very significant. All of the preceding results indicate 
that earthquakes and associated mass movements play a larger 
role in the excitation of the Chandler wobble than elastic 
dislocation theory (alone) predicts. More observations are 
needed before the earthquake excitation level can be 
accurately estimated. 

The earthquake related excitations may arise from: 

1. the main shock; 2. (pre- and post-seismic) movements of 
lithospheric slabs during subduction; 3. viscous relaxation 
(post-seismic) of the asthenosphere near the earthquake 
region; 4. (pre-seismic) lithospheric flexure (bowing) near 
subduction zones like that of Chile; 5. Aseismic fault slip. 
Subduction (zones where most earthquakes occur) appears to 
have a significant role in excitation of the Chandler wobble. 
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Chapter I . 


INTRODUCTION 

The Chandler wobble Is a polar motion of fourteen month 
period; specifically it is manifested as an angular displace- 
ment between an axis fixed to a stationary geometric pole on 
the surface of the earth and the spin axis of the earth. 

Since the discovery of the wobble in 189.1 (Chandler, 1891), 
the sources of its excitation have been extensively investi- 
gated without great success. J. Milne (1893) was one of 
the first to suspect some relation between polar motion and 
earthquake s. 

Extensive variation of latitude data have been recorded 
and analyzed to examine the driving source of the fourteen 
month wobble. More recently, theoretical elastic-dislocation 
models have been developed to determine if earthquakes can 
indeed drive Chandler wobble. Unfortunately, variation of 
latitude observations taken with zenith telescopes or other 
conventional optical techniques are not very precise; and 
e 1 as t i c -d i s 1 ocat i on models, while explaining immediate fault 
zone behavior fairly well, do not begin to model extensive mass 
movements associated with earthquakes or the inelastic 
behavior of the earth. This work examines the degree to which 
earthquakes and related mass movements excite the Chandler 
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wobble. 


Fortunately, we are not at an impasse in addressing this 
question. In the last two decades several techniques have 
been developed which more precisely measure the polar motion 
of the earth. These include Lunar Laser Ranging (LLR), 
Satellite Laser Ranging (SLR - via the LAGEOS satellite), 
and Very Long Baseline Interferometry (VLBI). VLBI is 
currently able to measure? the components of polar motion 
with an optimum precision of 0.5 mil 1 iarcseconds 
(Robertson, 1985). Since VLBI observing sessions span 
typically one to two days, one can readily study earth 
orientation on time scales up to these lengths. Using 
many day-long VLBI data sets, one can examine polar motion 
over time periods of weeks to years. This study looks at 
polar motion on time scales of several hours to several 
months . 

VLBI parameters such as polar motion components or 
baseline lengths (the chord distance between two points on 
the earth's surface) have been estimate d in the past from 
least squares analysis (Clark et al., 1985). A more power- 
ful technique for parameter estimation is Kalman filtering/ 
smoothing, and a major aspect of the work presented here is 
the application of Kalman filtering to VLBI observations to 
produce optimal estimates of VLBI parameters. 

The intent of this work is to investigate relationships 
between earthquakes (and associated mass movements) and the 
Chandler wobble of the earth using VLBI data. Some 
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questions to be addressed are: To what degree do earthquakes 
excite the Chandler wobble of the earth?; Does the Chandler 
wobble trigger earthquakes?; Do the movements of lithospheric 
slabs decoupled (from lithospheric plates) by earthquakes 
excite the Chandler wobble?; To what degree does unsteady 
bowing (lithospheric flexure) at ocean trenches excite the 
Chandler wobble of the earth? VLB I polar motion estimates 
appear to be approaching sufficient precision and temporal 
resolution to begin answering the preceding questions in some 
detail . 

The arrangement of the chapters in this work follows. 

In Chapter II polar motion fundamentals will be presented, 
while Kalman filtering and smoothing will be introduced in 
Chapter III. An introduction to Very Long Baseline Inter- 
ferometry will be given in Chapter IV, and a review of 
research into the relationships of earthquakes and Chandler 
wobble will be presented in Chapter V. In Chapter VI, a 
Kalman filter will be developed which produces polar motion 
parameter estimates from VLB1 group delay data. Filter 
optimization and operation will be explained in Chapter VII 
and VIII. Results and conclusions will be presented in 
Chapters IX and X. 
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Chapter II. 


THE EQUATIONS OP POLAR MOTION 

It is now appropriate to define some parameters of earth 
orientation and rotation. Consider that the earth is 
rotating with some angular velocity with respect to a quasi- 
inertial reference frame. The earth does not rotate at a 
constant angular velocity, and as such, there are changes 
in the length of day. In addition to rotational variations, 
there are small amplitude oscillations of the earth in 
directions orthogonal to the earth's mean rotation axis. The 
Chandler wobble makes up a part of these osc 1 1 1 at tons which 
are otherwise known collectively as polar motion. Chandler 
wobble differs from nutation because it is not a forced 
motion due to interactions of the earth with other planets, 
moons or stars (Stacey, 1977). 

The polar motion of the earth is made up of at least 
three motions: the Chandler wobble, the annual wobble and 
a secular trend in the polar motion. "The [Chandler] wobble 
results from rotation of the Earth about an axis that 
departs slightly from its axis of greatest moment of 
inertia." (Stacey, 1977). The Chandler wobble manifests 
itself as an oscillatory angular displacement between the spin 
axis and axis of greatest moment of inertia of the earth; 
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the period of oscillation is roughly fourteen months. 

Superimposed on the Chandler wobble is a wobble of a 
twelve-month period, the annual wobble; the annual wobble is 
driven by interactions of the oceans and atmosphere 
with the earth. The remaining component of polar motion 
discussed here is the secular trend. The secular trend is 
a long term drift in polar motion; it has a drift amplitude 
of several mi 1 1 iarcseconds per year (Morabito and Eubanks, 
1985). The secular trend appears to be driven by isostatic 
rebound of the earth after unloading by glacial melting 
(after ice ages). More details about polar motion can be 
found in Lambeck (1980). 

The changes in length of day alluded to previously can 
be measured by comparing time as kept using atomic time 
standards to time found astronomically (by looking at the 
rotational position of the earth with respect to quasars 
(these quasi- stellar objects define a quasl-inertial reference 
frame) or other stars). The latter time entity is known as 
sidereal time. The difference between atomic clock time and 
astronomical time is typically used to characterize changes 
in length of day. 

A type of astronomically-based time scale (like 
sidereal time) commonly used in VLBI work is UT1 (Universal 
Time). The atomic clock time scales, UTC (Coordinated Uni- 
versal Time) and TAI (International Atomic Time), are also 
used in the VLBI community. TAI differs from UTC by leap 
seconds which are introduced into UTC by the Bureau Inter- 
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nationale de l'Heure (B1H - a world time keeping agency). 

Leap seconds are corrections added into UTC to keep UTC 
within one second of time determined from earth rotation 
rate (Sullivan, 1984). The quantities UT1-UTC and UT1-TAI 
describe changes in the length of day. For further infor- 
mation about such time scales, please refer to Kaplan (1981) 
and Ma (1978) . 

The mathematical development which describes polar motion 
and rotation is now presented. Following the treatment (with 
modifications) of Nunk and MacDonald (1960), Conservation of 
Angular Momentum is used to describe the wobble and rotation 
of the earth 

dL 

N = + w X L . (11.1) 

dt 

N represents the torques acting on the earth, L is the total 
angular momentum of the earth, and u> is the angular velocity 
of the earth. We are concerned with the motion of the 
coordinate system xj (where i=l,2,3) which rotates with 
angular velocity w with respect to the inertial system Xj. 
xj and Xj are coincident at some reference epoch. 

It is now helpful to separate the total angular momentum 
L into two parts 

Li - C i j ( t ) oj j + 1 i ( t ) . (11.2) 

Cjj(t) is the inertia tensor and 1 ^ ( t ) represents a relative 
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angular momentum term. The components of C can be expressed 


as 


C 


i j 


P (x k x k 6 i j - x j x j ) d V 

JV 


(II .3) 


where p is the density of the object of concern, V is the 
volume, and 6jj is the Kronecker delta. The relative 
angular momentum components are given by 


li 


j P e ijk x j u k dV - 
V 


(II .4) 


where €ij k is the Levi-Civita symbol, and u k is the component 
of velocity with respect to the xj-system. 

Substitution of (II. 2) into (II. 1) gives the Liouville 
equat ion 


d 

N = (Cw+ 1) + u)X(Co)+ 1). (II. 5) 

dt 

Simplified versions of equation (II. 5) will now be formed 
to determine elementary equations of motion for Chandler 
wobble and changes in the length of day. 

One can also describe the Chandler wobble as a normal 
mode of the earth (Smith, 1977; Smith and Dahlen, 1981) and 
can use elastic-gravitational normal mode theory to charac- 
terize its behavior. It is an interesting note that 
Dahlen suggests "We could really say the Chandler wobble was 
just like every other normal mode if it were excited 
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principally by earthquakes, but that does not appear to be 
the case." 

Perturbation theory is applied to the Liouville 
equation (II. 5 ) to simplify it. For the case of polar motion 
and earth rotation in which movements of figure axes and 
rotation pole (and fluctuations in spin rate) are small 
with respect to the reference coordinate system Xj, the 
following perturbation scheme works well 


C 1 1 = A + C 1 1 

- c 22 

= A + c 2 2 . 

> c 33 = C + 

c 33 

c 12 = c 1 2 

- c 13 

= c 1 3 

■ ^ 2 3 = c 2 3 

(11.6) 

w 1 = Om 1 

, (jj 2 

= Ilm 2 

, U3 s (l+m3)0 

• 

In this model, 

A , A 

and C are 

the initial 

(pre-deformation) 


moments of inertia of the earth with C being the largest 
moment. Cjj are the perturbations to the inertia tensor. 
Clearly if all perturbations Cjj go to zero , we are left 
with an axi -symmetri c earth of major moment of inertia C and 
equal minor moments A. This makes sense and is a fair 
approximation to the geometry of the real earth. 

The angular velocities w = (u>j , W2> W3) are written as 
perturbations of the mean angular velocity of the earth, 0. 

and u >2 will be trivia] if m j and m2 go to zero, 
but cog is generally of order 0; this is intuitively 
gratifying, considering that one desires an earth which 
continues to rotate at some nearly constant angular velocity, 
regardless of the perturbations applied. 

Substituting definitions (II. 6) into the Liouville 
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equation (11.5), neglecting products and squares of pertur- 
bations and other small terms, one arrives at the following 
simplified Liouville equations 


m i m 2 

4 m 2 = 2 : - m i = - 0 1 (II. 7) 

Op O 


m 3 = 0 3 


( II .8) 


where 


Q^(C-A)® j = 0 2 c 13 + Qc 2 3 ■* 0.1 i + i 2 N 2 


0 2 ( C - A ) «l» 2 = n 2 c 2 3 " °c 13 + 0 1 2 - 4 Nj 


(11.9) 


o c $ 3 - *0 c 3 3 - o 1 3 4 n 


N 3 dt ' , 


where the dot indicates differentiation of the designated 
variable with respect to time, and o r = ((C - A)0)/A is a 
frequency parameter of the wobble. The r subscript denotes 
that a "rigid" earth is being modelled. Equations (II. 7) are 
very simplified descriptions of the Chandler wobble; they 
do not mode): the effects of an elastic earth on wobble; 
the lack of participation of the core in wobble; the effects 
of oceans; solid -fluid interactions at the core - mantle 
boundary, etc. For further explanations about the pre- 
ceding effects, please consult Lambeck (1980), Smith and 
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Dahlen (1981) and Dickman (1982). 

Equations (II. 7) govern the polar motion of our earth 
model, while equation (II. 8 ) characterizes earth rotation, 
mi , m 2 are the parameters typically measured and 
used to describe polar motion. The excitation functions 
<&1 , <&2 * are not simple, and their earthly sources 

are even more difficult to pin down. The atmosphere and 
earthquakes are possible sources of excitation, for example. 

The coordinate system first conventionally used for polar 
motion measurements was that defined by the International 
Latitude Service (II.S), an organization which made polar 
motion observations. The x 3 (polar) axis of the ILS 
system (determined from the mean latitudes of five ILS 
observatories from 1900-1905) is called the Conventional 
International Origin (CIO). The x i and x 3 axes are 
oriented along the Greenwich meridian and 90 degrees West 
meridian, respective ly, thereby completing the triad. 

By convention, m j is taken along the +- x j direction and m 3 
is in the - x 2 direction. 

This is a good opportunity to illustrate how the effect 
of an earthquake would appear in measured polar motion 
results. It is common to express the components of wobble 
(mj, m 2 ) as a complex number: m = mj * i m 3 

(also + 103 ). An equation for polar motion 

as a function of modified excitation function terms is 
presented in Munk and MacDonald (1960) 
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’1 

m(t) =■• exp(ia 0 t){m 0 io 0 4*( x ) exp ( i a 0 x ) dx } , (II. 10) 

- cw 

where mp is complex ( mp initial wobble amplitude), and 
Op is a wobble frequency term. 4* and op are modifi- 
cations of <K and o r , which account for the fact that the 
elastic earth deforms in response to its rotation. Thus the 
complex equation (II. 7) for a rigid earth, m - i a r ( m - ® ) 
can be extended to the case of an elastic earth 
m - iop(m - 4*) , where m, 0 and 4* are complex (Munk and Mac 
Donald, 1960; Lambeck, 1980). 

One can model the earthquake excitation 4> as a Heaviside 
step function, J H (t), 

0 , t<0 

JH( t) = < (11.11) 

J , t>0 

i 

where J is a complex constant. The step function is useful 
as an approximation of the earthquake excitation function 
because mass movements associated with the release of 
elastic energy during a seismic event should occur in a 
time short compared with the period of the Chandler wobble 
(Mansinha and Smylie, 1970). The earthquake occurs at time 
t^O. Thus for t < 0 

m(t) a mo exp(ioQt) (11.12) 

and for t > 0 
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jn(t)-m 0 exp(ioQt) + J(1 oxp(iaQt)) , 


( 11 . 13 ) 


There is no change in rotation pole position at t^O since 
m(t) is continuous. However, the pole about which the polar 
motion is taking place (the secular pole) does change as 
does the radius of the polar motion about this new pole, 
and thus m has a change in curvature. 

The effect of a hypothetical earthquake and its 
associated mass shift on polar motion is shown in Figure 1, 
taken from Mansi nha and Smylie (1970). A kink occurs in the 
path of the pole at the time of the earthquake, and the 
radius of the polar motion about the secular pole may also 
change. The kink is not a step discontinuity in the polar 
motion, even though the excitation function is discontinuous. 
It has been difficult to resolve kinks in existing polar 
motion data (Haubrich, 1970). It may be possible with Kalman 
filtered VLBI data in coincidence with a sufficiently large 
earthquake. The application of Kalman filtered VLBI data to 
the study of earthquake excitation of Chandler wobble is the 
central theme of this work. 

The Chandler wobble is driven by some mechanism (or 
mechanisms) and one would expect such an oscillation to be 
damped by other physical processes. Likely damping processes 
include mantle anelasticity and ”... a non -equil ibrium 
oceanic response to wobble.” (Dickman, 1986; see also Smith 
and Dahlen, 1981). The Chandler wobble quality factor, Q, 
which is used to characterize Chandler wobble (fractional) 
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energy dissipation, has been estimated by various techniques 
to be in a (preferred) range of 30 to 170 (Oickman, 1986). 

The value of Q is directly proportional to how long it 
would take to damp out the Chandler wobble. A Q of 100 
corresponds to a damping time of roughly 38 years. The 
excitations applied to the earth must be sufficient to 
drive the Chandler wobble to observed excitation levels, 
and one must account for losses in the driven earth system, 
due to damping mechanisms, when one is trying to evaluate 
the excitation sources. 
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Chapter III. 


KALMAN FILTERING AND SMOOTHING 

This section describes the merits and structures of 
Kalman filters and smoothers. Since the time of Gauss, 
scientists and engineers have used least squares analysis of 
measurements to provide best estimates of parameters of 
interest. Least squares is applicable when the statistics 
of an observable of interest are stationary. Data are 
handled en masse and a single set of parameter estimates and 
their statistical uncertainties is produced by the least 
squares analysis. This procedure works wonderfully for 
instructional phy s i cs- 1 abor a tory excercises in which one 
measures the dimensions of say, a metal block; but it leaves 
something to be desired in cases where dimensions are 
changing in time, such as finding optimal estimates of the 
distance and velocity of a rocket, which has just been 
launched. The Kalman filter algorithm was formulated in 
early 1960's (Kalman, 1960; Kalman arid Bucy, 1961) to deal 
with such situations. 

A Kalman filter is a recursive optimal estimator for 
parameters of interest. Rather than forcing statistical 
moments of observables to be stationary, the filter uses the 
modelled physical dynamics of a given problem to provide 
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optima] estimates of parameters of concern in the time domain. 
However, how optimal a parameter estimate turns out is 
dependent on the accuracy and completeness of the physical 
models employed in the filter. In general, measurements are 
handled sequentially, and from a series of measurements, the 
Kalman filter gives a series of parameter estimates in time. 
For instance, in tracking a rocket each range measurement 
input into a Kalman filter will produce an optimal estimate 
of distance and velocity of the rocket, from the filter 
shortly thereafter. This process can continue in a 
recursive manner until the payload has returned to earth. 

It should be noted that Kalman filters are limited in 
application to systems which satisfy the assumptions on 
which such a filter is based. The main assumption is that 
the Kalman filter be applied to systems which can be 
modelled as (mathematically) Markov processes (for details 
see Jazwinski (1970)). In addition, Kalman filters are 
particularly useful in describing stochastic processes 
such as white noise, random walks and integrated random 
walks. The filter may fail when other stochastic behavior 
operates; in such cases other optimal estimation techniques 
should be sought. 

The criteria for how optimal a Kalman filter is, are now 
defined. The optimal estimates of the filter or smoother 
should be unbiased, of minimum variance and consistent; that 
is, the estimates should achieve the true value of a 
parameter as the number of measurements increases (Gelb, 
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1974). The meanings of unbiased and minimum variance are 
evident . 

Before too much confusion is introduced by use of the 
words filter and smoother, let each be defined. A filter 
can produce estimates actively during the time over which 
data are being taken. A smoother is generally applied to 
data after the data have already been taken. Pictorial 
definitions of filtering, smoothing and prediction are 
displayed in Figure 2 (Geib, 3974). Both Kalman filters 
and smoothers will be developed in this work. 

It is instructive to begin by looking at the more 
familiar least squares analysis. An attempt will be made to 
follow the mathematical notation of Jazwinski (1970). In 
this dissertation (especially for Kalman filter notation) 
capital Roman letters represent matrices and small Roman 
letters are used to denote vectors. One can describe the 
linear measurement process of an experiment by the 
following equation 

y = M x t v (III.l) 

where y, x and v are column vectors and M is a matrix. In 
the equation above y is a (m by one) vector of measurements 
made in an experiment, x is a (n by one) vector of the 
parameters to be estimated (the ultimate results of an 
experiment) and M is a matrix (m by n) of factors (partial 
derivatives if the problem is non-linear) relating the 
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DENOTES SPAN OF AVAILABLE 
MEASUREMENT DATA 



(o) FILTERING 



(c) PREDICTION 


Fig. 2. Three Types of Estimation Problems 
(estimate desired at time t). (From Gelb et a 1 . 
1974 ) . 




measurements to the parameters, v is a (m by one) vector of 
measurement "noise" or uncertainty. n and ■ are positive 
integers . 

An optima], minimum variance estimate of x, namely x, is 
desired. Using the method of Gelb et a 1 . (1974), one finds 

x by calculating the value of x that minimizes the scalar 
cost function (net squared estimation error), J 

J s (y Mx ) T ( y Mx ) (III. 2) 

where the raised-T denotes matrix transpose. The result of 
this procedure is 

x - (MWy . (Ill .3) 

The -1 exponent refers to inversion of the matrix in paren- 
theses . Thus, the sum of the squares of the (y - Mx) 
elements has been minimized in accordance with the defini- 
tion that the optimal estimate be of minimum variance. 

Equation (III. 3) gives the best estimate of x for the case 
of unweighted least squares only. 

Equation (lll.l) is very close in form to the measure- 
ment equation of the Kalman filter. The discrete form of 
the Kalman filter is taken (as modified) from Jazwinski (1970). 
More complete descriptions of the filter can be found in 
Jazwinski (1970), Gelb et al. (1974) and Brown (1983). The 
equations which comprise the Kalman filter algorithm are 
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derived in Appendix I . 


The Kalman filter measurement 


equation is 

y k =M(k)x k +v k (III. 4) 

where y k , M(k), x k and v k are the same as for the least squares 
variables, except that the indices k refer to the time at 
which the observation was made. Consequently, t i me ( k *- 1 ) = t k f j 
is the observation epoch directly and temporally after 
time(k)~t k . For Kalman filtering, the time interval 
^k+1 ” t k does not have to be equal to the previous interval 
t k - t k l . All (III. 4) describes is the measurement 
of data vector y k , which is related to the parameter vector 
x k , with the noise vector v k influencing the measure- 
ment process at the observation time t. k . 

In contrast to least squares analysis, an assumption of 
stationarity is not required in the Kalman filter. Rather, 
the dynamical variations in the parameter x k are 
modeled using the state vet: tor model ( Jazwinski , 1970) 

x k + 1 = *(k+l,k)x k + F(k)w k (III. 5) 

x k+ j is the parameter at time t k + 1 and x k is the 
parameter at time t k . 0(k+l,k) is an (n by n) state 
transition matrix describing the deterministic variation 
in x going from time t k to t k+1 . w k is the state vector 
which models stochastic variations in x, and V ( k ) 
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is the matrix which relates w k to x. Equation (III. 5), 
as posed, thus models both the deterministic and stochastic 
dynamic variations of x in time. 

For completeness, the properties of the "noises" 
v k , w k are defined to be as follows (Jazwinski, 1970), 
where E{*} is the expectation value of * 

E(w k } = 0 , E { v k } = 0 ; 

E{w k w 1 T )=Q(k)6 kl , E{v kVl T }=H(k)6 kl (III. 6) 

and R ( k ) >0 . 

Again 6 k ] is a Kronecker delta, Q(k) is a system noise 
covariance and R(k) is a measurement noise covariance. 

The first four relations suggest that w k and v k are 
zero-mean, non-stat ionary white noise processes. The fact 
that R(k)>0 helps to assure that the Kalman filter does not 
diverge as it is applied. Equations (III. 4), (III. 5) and 
(III. 6) are the system model of the discrete Kalman filter. 

x(k]k) and P ( k | k ) are the Kalman filter analogs to 
average values and variances of parameters in least squares 
analysis. K ( k + 1 ) is the gain of a Kalman filter and is 
quite similar to the gain of an electronic amplifier. The 
notation (k|k), (k+l|k) which follows many of the Kalman 

filter symbols merits explanation. The (k|k) notation is 
simply that of conditional probability. x(k|k) is the 
optimal estimate of x at time t k conditioned on the 
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information in the previous measurements through yj< . 
x(k+l|k) is an optimal estimate at time t^j of the 
state parameter x based on the information of measurements 
y^; since measurement y^ + i is not included in the 
measurement y^. x(k+l|k) amounts to a prediction of 
of the estimate x at time t^i based only on measure- 
ments through y=y k . The notation is similar for error 
covariance P. Also, a small k in single parentheses 
following a variable, (k), denotes that the Kalman variable 
is at time t^. 

The operational steps for a discrete Kalman filter are 
(Jazwinski, 1970) 

(1) Store the state of the filter [ x ( k | k ) , P ( k | k ) ] ; (III. 7) 
where P(k|k) is the error covariance of x(k|k). 

(2) Project the filter state forward in time (a prediction 
step ) 

Sc ( k + 1 | k ) = 0 ( k+ 1 , k ) x ( k | k ) (III. 8) 

(3) Project the error covariance matrix forward in time 
(another prediction step) 

P(k + 1 |k)=<D(k+l ,k)P(k|k)® T (k + l ,k)+r(k)Q(k+l ) I’ T ( k ) (III .9) 

(4) Determine the Kalman filter gain matrix, K 
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( I [ I . 10) 


K ( k 4- 1 ) = P ( k + 1 | k ) M T ( k i 1 ) { M ( k + 1 ) P ( k + l | k ) M T ( k * 1 ) 

+ R ( k + 1 ) } ' 1 

(5) Combine measurement information y^fi with the projected 
state x ( k + 1 | k ) to form the optimal state estimate at 
time 

k(k+l | k+1 )«* (k+1 | k) +K (k+1 ) {y k+1 -M(k+i )S(k+3 | k ) > (III.lt) 

(6) Update the error covariance matrix 

P(k+l|k+l)*{I-K(k+l)M(k+l)}P(k+l|k){I-K(k+l)M(kM)}^(III.12) 

+ K ( k i 1 ) R ( k + 1 ) K T ( k t- 1 ) 

(7) Go to the next time step by setting k=k+l, and return to 

step (1) of this procedure. (III. 13) 

Filter steps (1) (7) are now explained. x(k|k) and P(k|k) 

are the most recent estimates of the state vector and its error 
covariance. They are stored in step (1). If one is just 
starting the Kalman filter x(k|k)=xg and P(k|k)=Pg, where Xg 
is the vector of initial state conditions and Pg is the 
initial condition of the error covariance. Usually xg and 
Pg can be crude guesses of x and P or they can be based on 
a priori information. If one has no a priori knowledge of 
the value of parameter x, one should specify Pg to be very 
large as an approximation to infinity. In filter steps (2) 
and (3), one uses the deterministic state transition matrix 
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<t> ( k <- 1 , k ) and the stochastic system (process) noise covari- 
ance Q ( k + 1 ) to predict (project forward in time) the state 
of the system in the future, { x ( k+ 1 | k ) , P ( k+ 1 | k ) } . This 
procedure relies on accurate system modelling (equation 
III. 5) to achieve consistent estimates. The filter utilizes 
all the physical information about the system to estimate 
state values and their uncertainty in the near future, i.e. 
the next time epoch tj< + j. <D = I (I is an identity matrix) 
means that x^ + i - x^ is constant in time (if no 
random forcing is assumed), and Q = 0 implies that the 
parameter vector xj< being estimated is strictly deter- 
ministic. 

In step (4) the gain matrix is computed. Clearly the 
gain combines the predicted covariance of the estimates 
P ( k -* 1 | k ) and the measurement noise covariance R(k + 1). The 
gain is the multiplier which determines how much the 
predicted parameter estimate x(k+l|k) will be adjusted when 
the next parameter estimate x ( k + 1 | k -* 1 ) is made. 

Finally, the updated parameter estimate x ( k •* 1 | k + 1 ) and 
covariance P ( k +• 1 | k <• 1 ) are calculated in steps (5) and (6) of 
the filter. These are the best estimates of x and P for time 
step tk + i- In step (5), [y^i - M(k+l)x(k^l|k)J is a quantity 

called the "innovation" which reflects the difference between 
measured data y^ + t and the modeled prediction of what the 
observation should be. Compare the innovation with the 
measurement equation (III. 4). Step (5) modifies the predicted 
parameter estimate x(k+l|k) by adding to it the product of 
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the gain and the innovation. Thus, the final state estimate 
at some time includes the most recent measurement information. 
Filter step (6) and its meaning are more complex than step (5). 
If the gain were zero, P(k+l|k+l) would simply be equal to its 
predicted value P(k+l|k). At the end of the filtering 
process, k is set to k+1 and the recursive filter loop 
begins again, and continues until the last measurement is 
processed. 

It is apparent that a matrix inversion is required in 
step (4) of the filter, and one would like to avoid such 
inversions if at all possible. Inversions take time to do 
and may fail computationally. In the case of application 
of the Kalman filter to VLB! delay data, matrix inversion 
can be avoided if VLBI delays are the only measurements 
utilized in y. The concept of a VLBI delay will be 
described shortly. It is sufficient to know that the delay 
is indeed the primary observable in VLBI data acquisition. 

Thus, for VLBI delay analysis, y^ becomes a scalar (m=l), 
and the matrix inversion of filter step (4) becomes a simple 
division. 

However, multiple delays are usually measured at a given 
time tfc , presenting another barrier to avoiding the inversion 
(m / 1). An "iterative" approach for handling multiple 
measurements is cryptically outlined in Jazwinski (1970). 

To implement this, one must first model the measurement 
noises v^ as uncorrelated. This allows one to process 
each delay at tj< in a one-at- a-time fashion through the 
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filter loop. Then, one proceeds through steps (1) (7) for a 

single delay, and thereby y^ is now most definitely a scalar. 
The procedure for avoiding the inversion follows. 

During the first run at time t^ through steps (l)-(7), 
the state transition matrix 0 ( k + 1 , k ) and system noise 
covariance Q(k+1) use their values as normally calculated 
(as calculated using the non iterative Kalman filter). 
However, upon completing this iteration, the time step is 
not increased from k to k+1. Rather, the filter is forced 
to stay at tinve t^. Now the second measurement at time t ^ 
is filtered (steps ( 1 ) - ( 6 ) ) , but 0(k+l,k) is set to 
identity (I) and Q(k) is given zero value. This is because, 
while the filter is iterating within multiple delay data at 
time tfc , the state vector x should not change determinist- 
ically (thus 0=1) nor stochastically (thereby Q = 0). 

Once all the delays at time t^ have been filtered, the logic 
of this paragraph is repeated at time t^ +1 , and so on. The 
iterative filter can thereby loop through uncorrelated 
measurement data without the added need of calling an inver- 
sion subprogram into use. 

The linear Kalman filter outlined in equations (111.7)- 
(III. 13) works well for linear physical problems. The 
procedure which makes the Kalman filter suitable for non- 
linear applications is termed M linearization 1 ' and while 
many authors do not describe it in the; literature because 
linearization supposedly is widely known, it will be 
outlined briefly here for completeness. 
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Continuous (as opposed to discrete) non-linear system 
equations analogous to the linear system equations (III. 4) 
and (III. 5) from Jazwinski (1970) are 

y(t k ) = h(x(t k );t k ) + v k (III. 14) 

dx £ 

and = f ( x t , t ) + G(t)w t . (111.15) 

dt 

y(t^) is the continuous measurement vector, h(x(t.^),tj c ) is a 
non-linear measurement function and is the measurement 
noise. x t is the continuous parameter vector, f(x t ,t) is the 
non-linear state dynamics function and Wf- is a white Gaussian 
noise of zero mean. G(t) is the matrix which relates the white 
noise vector to the parameter vector. 

Equation (III. 14) can be linearized by proposing the 
existence of a reference measurement y(t^) = h(x(tj < ),tj < ), 

where the bars over y and x denote reference values. In 
essence, one assumes reference values y and x which approx- 
imate the non-linear behavior of y and x. Subtracting the 
reference value y from both sides of equation (III. 14), one 
finds 

y ( t k ) -y ( t k ) =h(x( t k ) , t k ) -h( x( t k ) , t k ) +v k ( 1 1 1 . 16 ) 

If one defines the variation in the measurement y(t^) from 
the reference value y(t^) as 
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6y ( t k ) s y ( t k ) y ( t k ) 


( i i i . 17 ) 


and then assumes that the difference between actual and 
reference measurement functions can be represented by the 
first term of a Taylor series expansion 

8h i (x(t k ) , t k ) 

h(x(t k ) ;t k )-h(x(t k ) ,t k )* 6x(t k ) (111.18) 

8xj 


where 6x^_ s x^ x(t), 


( 1 1 1 . 19 ) 


then the new measurement perturbation equation is 

6y(t k ) = M{ t k ; x( t k ) }6x( t k ) + v k , (111.20) 

8hj (x( t. k ) ; t k ) 

and M { t k ; x ( t k ) } £ . (III. 21) 

8x j 

The variation in measurements 6 y ( t k ) is related to the varia- 
tions in parameters 6x(t k ) by the matrix of partials M. In 
awesome fashion, the non-linear equation (111.14), which 
describes behavior of the total measurement y(t k ), has been 
converted to an equation of Kalman filter measurement form 
(see equation III. 4). Equation (111.20) deals not with 
the total measurement y(t k ), but with the variation in the 
measurement 6y(t k ). Also equation (III. 20) is already 
in discrete form. 

Linearization of equation (111.15) is somewhat more 
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difficult. The reference value in x is given as 


dx ( t ) 


dt 


= f ( x ( t ) , t. ) . 


Subtracting (III. 22) from (III. 15) one finds 


( 1 1 1 . 22 ) 


dt 


(x t - x ( t ) ) = f ( x t , t ) f (x t , t ) +0 ( t ) w t . 


(Ill .23) 


Employing definition (III. 19), the relation 


af i ( x ( t ) . t ) 

f (x t , t ) -f (x ( t ) , t )* 6x t 


dx 


j 


( 1 1 I . 24 ) 


and defining 


F { t ; x ( t ) ) s 


9f t (x(t),t) 

3x< 


(111.25) 


(Jazwinski, 1970), one finds 


d(6x t ) 


dt 


= F { t ; x ( t ) ) 6x t + (i(t)Wf 


(111.26) 


(111.26) is still a continuous state equation and must be 
converted to discrete form. Details of the discretization 
procedure occur later in this work and will not be given here 
The discrete version of (III. 26) is ( Jazwinski , 1970 ) 


5x(t k + 1 ) = ®(t k + 1 ,t k ;x(t k ))6x(t k ) + w(t k + 1 ) (III. 27) 
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Equation (111.27) bears no small resemblance to equation 
(III. 5), but now the x variation, fix, is being modeled as 
opposed to x-total. 

The state equations (III. 20) and (111.27) are very much 
like (III. 4) and (III. 5). It can be shown (Jazwinski, 1970) 
that 


x(t k | t k )=x(t k ) + 6*(t k |t k ) (111.28) 

and that the covariances of x and 6 x are identical. The 
preceding facts allow one to use the linear Kalman filter 
(steps ( 1 ) - ( 7 ) ) to optimally estimate parameter variations 
and the covariance of these variations for the system of 
linearized equations (III. 20) and (III. 27). This lineari- 
zation scheme will work as long as the approximations in 
(III. 18) and (III. 24) hold true. One may have difficulty 
applying steps ( 1 ) - ( 7 ) to a strongly non-linear physical 
system, and other more sophisticated forms of filters such 
as the extended Kalman filter and adaptive Kalman filter may 
be needed (Jazwinski, 1970). Fortunately, the Very Long 
Baseline Interferometry (VLBI) state and measurement system 
equations are only weakly non-linear (Ryan, 1984) and steps 
( 1 ) - ( 7 ) are very applicable. 

The last segment of the Kalman filtering and smoothing 
chapter is devoted to the latter topic. VLBI parameters are 
not estimated in real time, and data from VLBI observing 
sessions are available for analysis only long after a 
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session is concluded. Thus it would make more sense to 


smooth VLB I delay data instead of only filtering it. 

The " fixed- interval " smoother is the Kalman smoother of 
interest here. It is applied to a total data span of 
fixed duration, typically one day long in the case of VLBJ 
observations. The time interval between data points is not 
necessarily constant. The following development of the 
smoother is adapted from Brown (1983). This f ixed - i n te r va 1 
smoother was developed by Rauch, Tung and Striebel (1965). 

Consider a data set containing N+l points. One applies 
the Kalman filter (steps (l)-(7)) in a forward direction 
(time progressing in a natural sense). As the process 
continues, the optimal filter estimates x(k|k), x(k+l|k) and 
the associated error covariances P(k|k), P ( k + 1 | k ) are stored 
in memory locations for use in the smoother. The need to 
store the estimates and covariances may limit the applica- 
bility of the smoother to problems whose x and P values can 
be accomodated in available computer memory. As the last 
datum is handled by the Kalman filter, the fixed-interval 
smoother is turned on. The smoother traverses the data set 
in a backward time progression. The smoother then smoothes 
the stored optimal estimates and error covariances. The 
filtering step is a forward pass through the data, and the 
smoothing step is a backward pass through the data. 

The first smoothing step uses the last estimates from the 
filter, (x(N|N), P(N|N)], as initial conditions in the 
following algorithm. The smoother gain A(k) is first calcu- 
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1 ated 


A ( k ) = P(k|k)* T (k+l,k)P 1 (k+1 |k) 


( I [ I . 29 ) 


Next the smoother opt imal estimate and error covariance 
are calculated 


X(kjN) = x(k|k) + A(k)<x(kH|N)-k(k+l|k)} (111.30) 

P ( k | N ) = P ( k j k ) + A(k) (P(k + 1 | N) -P( k + 1 | k) }A T (k) (III. 31) 

where k = N-l. N-2 ,0. Equations (II1.29)-(1I1.31) 

are the complete, discrete, f i xed - i n te r va 1 smoother algo- 
rithm. The smoother has no prediction steps like the Kalman 
filter; all the data which are to be processed are left from 
the filter. Another point is that the optimal estimates 
ft(k|N), P(k|N) are undeniably the result of smoothing since 
estimates at time t^ are based on N observations, as is 
indicated by the algorithm. Also note that a matrix inver- 
sion is present in the gain calculation, and that it cannot 
be avoided by iterative handling of the measurements as was 
possible in filtering. The state transition matrix 0 should 
be calculated as was done for the iterative Kalman filter, 
but note that the smoother traverses the filter estimates in 
reverse time order. 

The smoother state parameter estimate x(k|N) is some 
modification of the filter estimates x(k|k) by a factor 
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employing the product of smoother gain A ( k ) with some 
differential variable, which is the difference of a filter 
estimate and a smoother estimate, [x(k->l|N) ~ft(k + l|k)]. 

The smoother covariance P(k|N) is calculated similarly to 
x(k|N), with slight additional complexity because P(k|N) is 
a matrix, not a vector like x(k|N). While the smoother 
involves fewer calculational steps than the filter, prac- 
tical difficulties arise because one must store filter 
information for smoothing, and one cannot avoid performing 
matrix inversions while recursively smoothing. 

The smoothing estimates [ $c ( k | N ) , P ( k | N ) ] should 
generally be more precise than the filter estimates, 

[5t(k|k), P ( k | k ) ] . Figure 3 illustrates how the mean square 
error estimates of a Kalman filter and smoother may compare. 
Smoothing results in a smaller mean square estimation error 
than either forward or backward filtering alone, but as 
always there is a concomitant loss of resolution. In this 
work, the forward filter used with the VLBI group delay data 
will be iterative and linearized as has been described 
herein. The smoother will also deal with linearized varia- 
bles and will be of the f ixed- interval type producing opti- 
mal estimates of parameters of geophysical interest. 
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Fig. 3. Advantage of Performing Optimal 
Smoothing (From Gelb et al . , 1974). The 

abscissa is time (progressing) during the 
measurement process. T is the total 
measurement interval (in time). The 
smoot hing estimates are the result of optimal 
combination of the forward and backward filter 
estimates; this fact explains why the mean 
square estimation error increases during 
smoothing for times at the end of the 
m e a s u r e in e n t interval. 



Chapter IV. 


FUNDAMENTALS OF VERY LONG BASELINE INTERFEROMETRY 

Very Long Baseline Interferometry (VLBI) was a child of 
the 1960's, matured in the 1970's and is producing results 
in the 1980's. It is a radio astronomy technique which 
provides very precise values of parameters of astronomical, 
geodetic and geophysical interest. Using VLBI, distances 
between two points on the surface of the Earth can be 
measured to a precision of several centimeters over a length 
of thousands of kilometers. Positions of quasi-stel 1 ar 
objects (quasars) in space can be determined to about one 
mas, which is a factor of 10 to 100 better than other current 
astronomical techniques. With VLBI, one can look at contem- 
porary plate tectonic motions instead of several- mi 1 1 i on - 
year averages. Earth orientation angles and rotation rate 
can also be found to high precision. This chapter describes 
VLBI in some detail, so that one can understand the technique 
and its associated complexities. Many other reviews on 
VLBI are available in the theses and papers generated by the 
East Coast VLBI Group (U.S.A). For example, see Robertson 
(1975), Ma (1978), Herring (1983) and Lundqvist (1984). 

Much of the VLBI work which is germane to geodesy, plate 
tectonics, and tectonic plate stability is performed under 
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the auspices of the NASA Crustal Dynamics Project (CDP). 

VLB I measurements will be regularly made via the CDP until 
at least the year 1990. VI,BI is occurring between points 
on the North American continent, Hawaii, Japan, Europe, 
Kwajalein, etc. The majority of these observations are 
carried out by groups affiliated with the East Coast (U.S.A.) 
VLB I group or NGS (National Geodetic Survey). 

Earth orientation VLBI measurements are carried out 
at least once every five days by the IRIS (International 
Radio Interferometric Surveying) network with stations in 
Massachusetts, Florida, Texas, Sweden and Germany. The 
measurements of the IRIS network have effectively been 
carried out since early 1984. The IRIS data set of 1984- 
85 will be the primary object of analysis of this research 
work . 

Simply put, VLBI is an interferometric experiment over 
a very long baseline (or baselines), which is usually a line 
between two surficial points on Earth of length up to about 
10,000 kilometers. The signals for the interferometry 
are radio wavelength and are generated by quasars and other 
distant radio sources. A radio signal arriving at the earth 
from such a source can be well approximated as a plane wave. 
The radio signals can be received by radio antennas and 
recorded with precise time information. If the respective 
arrival times at several antennas of a given wavefront from 
a radio source can be determined, then one can find the 
difference in arrival time at one station with respect to 
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a ri o t. her. 


This difference in a r r i v a J time is called the delay, and 
results from one VLBI antenna being closer to the source of 
the plane wave than the other antenna. The delay provides 
information on the distance between the two stations, and on 
t h e orienta t i o n o f the b asnlino with respect to the source. 
"The technique can be likened to the manner in which two 
blind people far apart on a beach can learn the direction 
from which waves are coming. If they recorded the exact 
time when each wave reached the toes of each observer they 
could determine the angle of its arrival." (Sullivan, 1984). 
They might also begin to know the distance between themselves 
after the arrival of many waves. Of course the previous 
thoughts are severely simplified. 

The simple geometry of a single baseline delay measure- 
ment is shown in Figure 4 (Ma, 1978). The baseline B 
connects the two VLBI stations and S is the unit vector 
pointing from earth to the radio source. U is the distance 
most directly associated with the delay being measured; an 
incoming plane wave traverses distance D during the delay 
time. The plane wave arrives first at station one, which 
is the VLBI station closest to the source. The plane wave 
then arrives at station two. The radio signals are 
recorded independently on t: a p e recorders at each station. 

A clock, usually a hydrogen maser time and frequency 
standard, is used to determine the time of plane wave 
arrival at each station. The tapes from each station are 
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Fig. 4. Basic VLBI geometry (from Ma , 1978 ). 
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sent to a "correlator" where the delay between the arrival 
times is found precisely. VLBI differs from connected- 
element interferometry in that the endpoints of the inter- 
ferometer are not connected electrically. A maser at each 
VLBI station provides time information on the data tapes to 
permit comparison of signal arrivals long after an experi- 
ment is run. 

The geometrical delay, T g , associated with the 
station geometry (see Figure 4) is 

-B • S 

T g = (IV. 1) 

c 

where c is the velocity of the plane wave and the dot 
denotes a dot product of two vectors. The geometrical 
delay is simply the distance D travelled by the plane wave 
between the two stations divided by wave velocity. The 
angle between vectors B and S can be estimated using 
relation ( IV . 1 ) . 

While the fundamental physics of VI.BI is trivial, 
real world phenomena make VLBI modelling complex. The earth 
rotates as each delay measurement is made; so the VLBI 
observation geometry is not constant in time. The geometric 
delay changes temporally, and the motions of the stations 
cause Doppler shifts of the incoming radio signals. The 
earth wobbles with annual and f our teen- month periods, and 
precession and nutation occur. The earth is not rigid 
in general, and effects of earth and ocean tides on VLBI 
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observations must be modelled. The incoming radio waves 
traverse the ionosphere and troposphere, and propagation 
through these media affects the measured delays greatly. 

Most of these effects are well described in Ma (1978), 
Robertson (1975) and Lundquist (1984). 

Herring et al . (1985) give a very clear description of 

several coordinate systems used in VLBI work and how one 
transforms a vector (such as a VLB] baseline vector) between 
the coordinate systems. The reference time epoch used in 
VLBI research is based on the year 2000, and is known as 
J2000.0. The time epoch formerly used in VLBI work was 
1950. The notation of Herring et al. (1985) differs 
somewhat from that used in this dissertation. b j is 
B 2000 herein; b t (Herring) is B terrestr1al 
herein. The polar motion components Xp and yp are 
X BIH anc * ^BIH herein. The mean obliquity of data 
(f.0 Herring) is given herein as G. Little else is 
different. The passage from Herring et a.l . (1985) follows. 

In this appendix we discuss in detail the rota- 
tion matrices which are used in our analysis to 
transform between 'crust-fixed' and inertial coor- 
dinate systems. Our terrestrial coordinate system 
is defined by a set of site coordinates which would 
be invariant with time if there were no tectonic 
motions of the VLBI sites. From these site coordi- 
nates we compute the vectors between the sites, b c , 
in the crust-fixed coordinate system. For the central 
epoch of each observation (typical duration 100400 
seconds), we compute the displacements of the sites 
due to the solid earth tides using the algorithms 
discussed in Herring et al. 1981. These displace- 
ments are added to b c to yield b t = b c + 
u e , where u e is the (three dimensional) earth 
tide displacement computed from a planetary ephemeras 
of the sun and the moon, and the Love numbers h=0.609 
and 1=0.085. (We assume zero phase lag.) We now 
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apply a sequence of rotations to the coordinate system 
in which b t is given to determine the coordinates 
of the baseline in inertial space, bj. This 
sequence of rotations is: 


bj[ - PNSW b t . 


(A. 1 ) 


We now discuss each of these rotations. We give 
firstly the expression used in evaluating the elements 
of the rotation matrices. We then discuss the 
physical significance of the transformations. In 
our discussion we will use the standard rotation 
matrices (Goldstein 1950, p. 109): 



1 0 0 \ 


1 C 0 - s ' 

R X U) = 

0 c s 

, RyU) = 

0 1 0 


O S c 
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\ S 0 c j 


( c s o' 
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* 
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o 



where C 

= c o s £ , 

and S = s i n £ 

, and £ 


of the rotation. 

We discuss the rotations in the order- 
are performed. The W matrix is given by: 


the argument 
in which they 


W - R x ( -y p )Ky(x p ) , (A. 2) 

where x p and y p describe the pole position 

(see later discussion) with the sign convention used 

by the BIH. The S matrix is given by: 


S * R z ( --GAST ) , (A. 3) 

where GAST is Greenwich apparent sidereal time given 
by : 


GAST = GMST 0 + [ cl ( GMST ) / d ( UT ) ] UT1 + A4>cos G, (A. 4) 

where GMSTq is Greenwich mean sidereal time at 
0 hr UT, d ( GMST ) / d ( UT ) is the derivative of GMST with 
respect to universal time, UT (Aoki 1984), and A^cosG 
is the 'equation of the equinoxes' with A^ being the 
total nutation in longitude of date and G being the 
true obliquity of date. The N matrix is given by 


N = R x (-£ 0 ) R z( A^Rx^O^G) , (A. 5) 
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where 6 0 is the mean obliquity of date, and and At 
are the nutation angles computed from the IAU 1980 
nutation series (Wahr 1981, Seidelmann 1982). The P 
matrix is given by 

P = Rz (?, a ) Ry ( ~®a ) R z ( z a ) * (A. 6) 

where £, a , 0 a , and z a are the standard 

arguments for precession (Lieske et al. 1977). 

We now discuss the meaning of the transformations 
given above. The W rotation moves the pole T of the 
terrestrial coordinate system to the pole R of an 
intermediate coordinate system (see Figure A.l)[Fig. 

5 ]. The pole R is approximately coincident with the 
pole of the instantaneous rotation axis (we discuss 
later in the appendix the specific definition of R). 

The S matrix rotates the coordinate system about the 
R pole. The combined PN matrix moves the pole R to 
the pole I which would be fixed in inertial space if 
the arguments of all of the transformations were 
correct. If we use the IAU [International Astronomi- 
cal Union] 1980 nutation series in the evaluation of 
N, the PNSW transformation is not exactly correct 
because we perform the sidereal rotation about the R 
pole rather t h a n * a b o u t the earth' s^instantaneous 
rotation axis, R However, the R pole 
moves relative to the crust with approximately a 
diurnal period under the action of the external luni- 
solar torques and, hence, in addition to the W trans- 
formation (whose arguments are assumed to change 
slowly compared to the diurnal period), we would need 
to apply another transformation, the diurnal -po lar- 
motion (I)NP) transformation (Kinoshita et al . 197§). 

Through second order in the separation of R and R 
(maximum separation 25 mas), this latter technique 
yields the same transformation as the technique of 
rotating about R (see discussions in Kinoshita et al . 
1979 and Seidelmann 1982). Consequently, the axis for 
which the nutations are calculated is one which would 
be fixed relative to the crust of the earth if the 
earth were not undergoing free motion. This axis is 
the axis of figure of the Tisserand mean outer surface, 
or simply, the 'body axis' (Wahr 1981). The inter- 
pretation of Xp and yp is more complicated 
than that for the nutation angles because the polar 
motion is composed of two types of motion: the free 

motion (Chandler wobble), and the forced motion 
(which has an approximately annual period and is due 
to interactions between the solid earth and the 
atmosphere/ocean system). 

Herring et al. (1985) make corrections to the 1980 IAU 

Nutation Theory. A clear account of the physics of 
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Fig. 5. The geometry of the various poles used 
in the transformation from the 'crust-fixed' 
coordinate system to the inertial coordinate 
system. The figure is not to seaie. The sepa- 
rations of various poles are: T and R« 300 mas, 
t R and 1 « 9,000 mas (due to nutation) p .1 y s « 0.1 

| degree (due to precession), and R and R 

*25 mas. (From Herring et ai . , 1985). 
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precession and nutation can be found in Stacey (1977). 

Now that the various transformation matrices and 
reference frames have been formulated, one should have 
sufficient courage and background to be exposed to the 
actual delay expression, which differs from the simple geo- 
metrical delay T g (FV.l), due to Doppler and rotation 
effects etc. It. is (Robertson, 1975) 

X =t 0 R- {R 2 -Ri>- <R R 2 + R R2> t <W LliT ) x 0 (IV. 2) 

where 

T 0 - (Ri K 2 ) * S - { (R+R 2 ) ■ S) { (R r R 2 ) ' S ) { 1 - (R + R 2 ) * S} (IV. 3) 

-< 1/2) { (R*R 2 ) ’ S){ (Rj -R 2 )- S> 2 + t A . 

The variables of (IV. 2) and (IV. 3) are defined here (Ma, 

1978 ) 

x - the delay (usually in nanoseconds) 

R = solar system barycentric position of the earth's 
center 

R I = geocentric position of station one (see Figure 4) 

R 2 = geocentric position of station two 
S = unit vector directed from earth to source 
x A = propagation media delay 
LPT = = Long Period Terms (explained in Ma , 1978). 

Dots over variables indicate differentiation with respect to 
coordinate time. Dots between variables denote dot products. 
Expressions for the delay rate, t, are not quoted here because 
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they will not be used in this work; the delay rate equations 
can be found in Ma (1978). The value of the velocity of light 
is set to unity (c=l) in expressions (IV. 2) and (IV. 3). 

From (IV. 2) and (IV. 3) it is evident that the delay 
t is directly proportional to the simple geometrical delay 
(IV. 1), which is ( R j - R 2 ) * S , because baseline length B 
~ R 2 “ R 1 • Also the delay is directly proportional to 
the propagation media delays, which makes sense. 

Equations (IV. 2) and (IV. 3) include terms resembling 

’ * * • ( f t ♦ 4 

velocities (R,R 1 ,R 2 ) and accelerations ( R , R j , R 2 ) 

which are due to a Taylor expansion procedure in the 

derivation of the delay. For the details see Robertson 

(1975). Equations (IV. 2) and (IV. 3) are approximations 

- 1 5 

which neglect terms of order less than 10 seconds in 
delay (Robertson, 1975). This approximation is reasonable 
because group delay residual standard deviations are 
generally on the order of 100 picoseconds or less. 

Relatively little has been said about the general method 
of determining the observed delay. For the moment, let the 
group delay be designated by the symbol At. x will be used 

again later to represent the group delay. The delay between 

plane radio wave arrivals at two VLBI stations can simply be 

expressed as (Shapiro and Knight, 1970) 

cAt =n\+6\;0<6\<\. (IV. 4) 

where c is the velocity of the radio waves, n is an integer, 
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X. is the wavelength of the radio wave of concern, and 6X. is 
some fraction of a wavelength. Note that the case of a non- 
dispersive radio wave is being treated here. Clearly, the 
propagation distance, cAt , between two stations can be 
thought of as a combination of integral and fractional 
wavelengths, as (IV. 4) suggests. Using the relation 
f = c/\, where f is frequency, one can convert (IV. 4) to 

2trf At = 2 nil + 4> (IV. 5) 

where 4>, the phase, is given by 0 - ( 2 rc 6 X. ) / \ . 

Differentiating (IV. 5) with respect to frequency, one 
can readily find (Shapiro and Knight, 1970) 

d<D 

= 2 K At . ( I V . 6 ) 

df 

Thus if one has phase information as a function of radio 
frequency, one can readily determine the group delay, At. 

The slope of the phase vs. frequency determines the group 
delay. VLBI observations are commonly made at multiple 
frequencies within two wide bands of radio frequency known 
as (Lundquist, 1984) S-band (2.3 Giga Hertz) and X-band 
(8.4 GigaHertz). The multiple frequencies within each band 
are given in Lundqvist (1984). The technique of determining 
group delays from phase observations at multiple frequencies 
within wide radio bands is called "frequency synthesis”. 

The technique is due to Rogers (1970). The reason for 
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making VLBI observations in both the S and X bands will be 
explain e d subsequently. 

From (IV. 6) one can derive an expression for the error 
in the group delay, 6At (Shapiro and Knight, 1970) 

6<J> 

6 At * ( I V . 7 ) 

2* ( f max " * mi n ' 

6® is the error in the phases 0 and f max , f m j n are the 
maximum and minimum frequencies enclosing the frequency 
synthesis bandwidth (Shapiro and Knight, 1970). From 
(IV. 7) it is evident that the error in the delay, 6At , is 
independent of the baseline length between two VLBI 
antennas . 

There is considerable work done in preparing and 
executing a VLBI experiment, so it will be instructive to 
describe the process. The logistics of the experiment, such 
as equipment and personnel dispatching, etc., will not be 
mentioned. More than two stations usually participate in an 
experiment. The station antennas observe, in unison, about 
ten or more radio sources sequentially up to a total of 200 
scans in 24 hours. In NASA Crustal Dynamics Program experi- 
ments, observation schedules are designed using the GSFC 
program SKED with the guidance of a VLBI scientist, and these 
schedules are sent in advance to the participating VLBI 
stations. Schedules must advise each VLBI station when to 
start and stop viewing a source, and these observational time 
limits must reflect the times when a given source is visible 
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at each station. Sources must be scanned long enough to 
produce acceptable interference fringes. There are many 
other complexities not addressed here. 

The signals received at an antenna are stepped down from 
radio to 2 MHz baseband frequencies and the resulting signals 
are stored on tapes. Each tape also has time information from 
the station maser. Once the experiment is completed, all 
the tapes are sent to a correlator location. Crustal 
Dynamics Project tapes are sent to the Haystack (Massachu- 
setts) correlator. The correlation process and much of what 
else goes on in VI.BI experiments are described in Ma (1978). 
The observed delay for each baseline at each scan is the 
ultimate output of the correlator. 

The delays are then sent on tape to an analysis team, 
such as the VI.BI group at Goddard Space Flight Center (GSFC). 
Theoretical delays for each baseline are calculated using 
the program CALC. The theoretical delays are used as 
reference values in a linearization scheme similar to that of 
equation (III. 17) and are based on a priori station 
positions, quasar positions, UT 1/polar motion values, etc. 

The difference between the observed delay and the theoretical 
delay is called the delay residual. Typically group delay 
residuals are analyzed, but phase delay information is 
avai lable . 

The delay residuals are used to solve (either by least 
squares analysis or Kalman filtering/smoothing) for adjust- 
ments to the a priori parameter values of station positions. 
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earth orientation angles, clock variations etc. The esti 
mation program at GSFC is called SOLVE* and generally a 
SOLVE solution requires several hours or more of analyst 
time to process a single VLB! experiment. Much of the 
analyst’s time is used deciding polynomials to approximate 
the clock variations (arising from station maser behavior 
and other physical phenomena) seen in VLBI residual data. 

One goal of this work is to see if a Kalman filter can 
solve for clock variations in the time domain with limited 
analyst intervention, thus reducing analyst effort and also 
eliminating the subjectivity of analyst clock modelling. 

VI.B1 solution results generated by SOLVE are transferred 
into various data bases, one of which is openly available to 
geophysicists and astronomers. Outside users can use the 
information to study current plate tectonics, crustal block 
stability, earth rotation and orientation, or radio source 
structure and position. 

The conversion of delay residual information to VLBI 
parameters of physical interest generally progresses towards 
solutions with better fits. Many SOLVE analysis sessions 
are required before solutions achieve stable form. In 
early (Initial) solutions* B1H rapid service ( prel i minary ) 
UT1 and polar motion values are used as a priori reference 
parameters. Later, final values of UT1 and polar motion are 
taken from the BIH Circular D. Meteorological information 
for calibrating the atmospheric delays in the observed 
delays is Included in SOLVE runs as soon as it is available. 
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If the physical models in CALC or SOLVE are improved, new 
computer runs are made with updated versions of the 
software. The process continues until a solution fit cannot 
be improved, currently less than 50 picoseconds RMS group 
delay residual error. 

There are some physical phenomena which grossly affect 
delay data. The models used by the East Coast VLBI group to 
minimize the effects of these phenomena on delays are 
described below. These include: the ionospheric correction, 
"clock" models, cable calibration, and various atmospheric 
models. This section is intended as a brief and contemporary 
explanation; some other phenomena affecting delay values, 
such as antenna axis offset, are not covered in this work, 
and generally can be found in Ma (1978). 

It was mentioned previously that VLBI delay data are 
taken in two frequency bands, the S and X radio bands. This 
is done to remove the effects of the ionosphere on the delay. 
Charged particles are present in the ionosphere and delay 
incident VLBI signals. Delay due to passage of VLBI signals 
through the ionosphere (Tj on ) is frequency dependent 
and can be modelled via the relation (Ma, 1978) 

T ion » k/f 2 . (IV. 8) 

where k is a constant and f is the observing frequency of 
the band of interest. 

Following Ma (1978), if and x 2 are group delays 
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observed at frequencies f* and f 2 (for instance, at X and S 
band) and t actua j is the actual group delay with ionospheric 
delay removed, then one can write 


2 

x 1 “ x ac tua 1 h k / f l 


(IV. 9) 


x 2 


x actual 


+ k/f 2 2 


( IV . 10 ) 


Since tj, 1 2 , fj and f 2 are known for a VI.BJ experiment, while 
k and t act - ua j are unknown, there are two equations and two 
unknowns. One can solve for k and ^actual- With some 
slight algebraic manipulation 


Ul t 2 )(fj) 2 (f 2 ) 2 
(f 2 ) 2 (f x ) 2 


(IV. 11 ) 


and 

< X 1 - x 2 ) 

T actual = x l" 2 (IV. 12) 

1 - (fi/f 2 ) 

r ac tual is a single delay deduced from observations at two 
frequency bands and is relatively free from any ionospheric 
delay contribution. 

Once a radio signal has penetrated the earth's iono- 
sphere, it must traverse the atmosphere. As the VLBI signal 
passes through the atmosphere, the refractive index of the 
medium increases, and the path of the plane wave is no longer 
straight. Ideal and actual radio signal paths through the 
atmosphere are illustrated in Figure 6 ( Lundqvist , 1984 ; 
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Elgered, 1983). Clearly the actual curved path travelled 
by the radio wave is longer than the ideal straight path that 
would be taken if the earth had no atmosphere; thus there is 
a delay associated with the propagation of a VLBI signal 
through the atmosphere due to additional path length. The 
atmospheric delay is roughly due to a dry component, from 
the air in the atmosphere, and a wet component, from 
water vapor in the atmosphere (an approximate description). 

The additional path length due to the dry component, L d , 
is given by (Moran, 1976; Ma, 1978) 

R 

L d - 77.6 P 0 (cm), (IV. 13) 

mg 

where 

7 

R = 8.3144 x 10 erg/mole = universal gas constant 

2 

g = acceleration of gravity (at station) - cm/sec 

m ~ molecular weight of dry air - 28.966 

Pq ~ dry air pressure (millibars). 

It is much more difficult to estimate the wet component 
path length because of the highly inhomogeneous distribution 
of atmospheric water vapor. Consider that the determination 
of the moisture content in the atmosphere plays a large role 
in weather prediction, which is problematic at best. One 
might then appreciate the task of estimating the wet compo 
nent path length faced in VI.Bl. The wet component path 
length depends on the amount of water vapor in the direction 
of the radio source being observed. A means of determining 


53 


the wet delay observational ly using a Water Vapor Radiometer 
(WVK) will be set forth later. The wet component path length 
is roughly ten to fifteen percent of the dry component path 
length (Ma. 1978). 

Even if directly observed atmospheric delay information 
is available, one normally adjusts the atmosphere delay in 
the zenith direction at each VI.B1 station. A priori atmo- 
spheric delays can be determined using barometric pressure 
(air column height) calculations or they can be set at 
atmospheric delay values from previous VLSI experiments. 

It is desirable to determine atmospheric path delays 
independently of VLBI measurements by using WVR's or radio- 
sondes, but until recently this has rarely been done. 

Several scaling equations relating zenith delay to 
actual path delay are, or will be, commonly used by East 
Coast VLBI analysis teams. The first is due to C.C. Chao 
(1972), and is a modified cosecant function (Elgered and 
Lundqv i s t , 1984 ) 


^ z 

L(EI.) *= , (IV. 14) 

sin(EL) + 0.00143 

tan ( EL ) + 0.0445 

where L(EL) is the path length at the elevation angle EL and 
where L z is the zenith path length. Elevation angles are 
measured relative to the horizon, so that the EL of zenith 
equals 90 degrees. The Chao model constants were found by 
best least square fits to radiosonde profile ray-tracing 
results (Elgered and Lundqvist, 1984). The Chao model only 
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requires eievation angle and zenith delay length for appli- 
cation, so it is generally applied to VLBl observations in 
which no ancillary meteorological data were recorded. 

A more sophisticated mapping function is due to Marini 
(1974) for elevation angles greater than ten degrees (Elgered 
and Lundqvist, 1984) 

1 ( A + B ) ( cm ) 

L = (IV. 15) 

f ( 1 a t , h ) s i n ( E L ) + B 

(A + B) ( s i n ( EL ) +0 . 0 1 5 ) 

where 

A = (2.277x10~ 3 )x(PM 1255/T +• 0.05) x e) 

B = (2.644x10 3 )x exp ( -0 . 14372xh ) 

h - height of station above sea level (kilometers) 
lat = station latitude 

f ( 1 a t , h ) - 1 . 0 - (0.0026 x cos(2 lat)) - (0.00031 x h) 

e - water vapor partial pressure (millibars) 

P - total pressure (mbar) 

T =* absolute temperature (K) 

Obviously meteorological data must be recorded at the 
V L B [ observing station in order to apply the Marini model. 

If meteorological data are available, the Marini model is 
generally used in preference over the Chao model for atmos- 
pheric delay calibration. However, VI«BI analysts are free to 
choose the atmospheric model that allows the best least 
squares fit in the delay residuals. 

Davis et al. (1985) give a newer relation between zenith 
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path length and path length at a specific elevation 

L z 

L ( EL ) ( IV . 16 ) 

sin (EL) + a^ 

tan (EL) + a 2 

sin(EL) + 33 

where a^ , a2 . a 3 are again determined from best fits to ray 
tracings. This mapping function gives much better results 
than the Chao and Marini models, especially at elevation 
angles below about ten degrees (Davis et al . , 1985). 

As was suggested earlier, it would seem better to observe 
atmospheric delay directly, and remove it from VLBI data, 
than to find the zenith atmospheric delay from VLBI data. 
Since the dry part of the delay is readily determined, it 
is only necessary to measure the contribution of the wet 
atmospheric delay at elevation. The wet delay can be calcu- 
lated by inverting brightness temperatures observed at 22 
GHz and 31 GHz using water vapor radiometers (Clark et al . , 
1985; Resch et a 1 . , 1984). Brightness temperature is a 

quantity which characterizes the radiation given off by a 
black body using Planck’s radiation law; in this instance 
microwave radiometers scan the radiation given off by the 
atmosphere. The water vapor radiometer technique relies on 
the fact that a water vapor resonant absorption line exists 
at 22 GHz, but not at 31 GHz. Thus the integrated water 
vapor content along some elevation angle can be estimated 
from the brightness temperature results. Several water 
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vapor radiometers are being tested at Crustal Dynamics 
Project (CDP) VLBI stations and should see active duty in 
the next several years. 

Since WVR data have rarely been availabJe in the past, 
it would be useful if some means could be discovered to 
retrieve wet atmospheric delay results directly from VLB! 
observations, so that VLBI data taken in the past decade 
could yield better results. Using a Kalman filter, Davis 
et al. (1985) were able to estimate wet atmospheric delays as 
a function of time (at each VLBI observation). These atmos- 
pheric delays, formed by Kalman filtering VLBI delay data 
only, match well (aside from slight biases) WVR atmospheric 
delay results observed independently of the VLBI delays. 

VLBI delay observations also include the effects of 
"clock" variations. Such clock variations are time-like 
signals which are superimposed on VLBI delay signals; the 
clock variations are separated from VLBI signals through 
least squares or Kalman filter analysis. Differential 
clock variations between stations add directly into the 
observed VLBI delays, and are principally due to three 
sources: maser time (and frequency) variability, phase 
calibration and cable calibration signal failures. 

The maser at a VLBI station is used as a time and fre- 
quency standard. A precise frequency standard is needed to 
step the incoming radio signals down to video frequencies 

suitable for recording. The hydrogen masers are typically 

1 4 

stable to one part in 10 in time, with hoped-for 
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improvement by one or two orders of magnitude as supercooled 

maser cavities come into operation. However, masers behave 

less well if their environment is not thermally controlled 

or if they are not isolated from other ills such as stray 

magnetic fields. The nuisance signal put into a VLBI 

delay by two masers at each end of a given baseline is the 

relative time difference between the masers. Most of 

the time difference between masers is a "ramp" function in 

time (otherwise known as a "rate"). Such a constant change 

in maser time as a function of time is due to an offset in 

the frequency of each maser, and can also be observed in a 

controlled laboratory environment. More often than not, any 

non-linearities visible in a maser-maser intercomparison are 

due to electronics external to the masers; the masers are 

1 4 

indeed stable to one part in 30 in time (Chiu, 1984). 

Laboratory maser behavior indicates that "clock" fine 
structure irregularities usually seen imposed on the 
dominant clock rate trend are created by phase calibrator 
signal variations or other changes in signal path. The 
phase calibrator at a VLBI station generates pulses from 
the station frequency standard and injects them into the 
main signal path at the feed horn of the receiver. The 
phase calibrator pulses then follow the same signal path as 
the quasar signals until their recording on magnetic tape 
in video format. The purpose of the phase calibration is 
to correct for retardation and dispersion in the cable and 
electronics as the signal passes from feed to recorder. "The 
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arrival tine associated with a signal is derived directly 
from the phase calibrator injected with that signal [the 
signal from maser to phase calibrator!." (Ma, 1978) 

Unluckily, the electrical reference signal passing 
from the station frequency standard to the phase injector 
is also carried by cable, and is subject to environmental, 
dispersive and retardation effects, thus changing measured 
delay values. Care must be taken to maintain the integrity 
and physical dimensions of the maser-to-phase injector 
cable. Several members of the East Coast VLBI group are 
currently working to minimize phase calibrator (and asso- 
ciated) variability affecting the total VI.BI delays. Such 
variability ultimately shows up in VLBI delay residual data. 
Cable calibration techniques are used to control the effects 
of thermal expansion and cable stretching on the length of 
the cable running from the station frequency standard to the 
VLBI receiver(s). 

Together, maser time fluctuations, phase and cable 
calibration related alterations, and other jumps are 
referred to as "clock" variations by the VLBI analyst. 

Clock variations must be removed from VLBI delay data in 
order to obtain reliable astronomical and geophysical 
results. This is done within the VLBI program SOLVE with 
considerable analyst intervention. The analyst looks at 
delay residuals plotted as a function of time. (Note that 
the process being described here takes place only for least 
squares analysis, not in Kalman f i 1 ter ing/smoothing. ) 
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Trends in the residuals are usually apparent, and the 
analyst tries to fit the delay residual plot with a poly- 
nomial (or polynomials) of the form 

T c = a 0 + a l* + a 2* 2 + ••• + A sin(ftt) + B cos(ftt) (IV. 17) 

where T c is the clock polynomial, a 0 . aj, a 2 A, B are 

constants which are evaluated in SOLVE to provide the best 
least-squares fit to the delay residuals, t is some measure 
of the time during the observing session, and ft is the 
diurnal rotation frequency of the earth. Sinusoidal 
(diurnal) behavior Is sometimes observed in the residual 
plots, and this can be fitted using the constants A and B. 

Sometimes one cannot fit a single polynomial to the 
clock variability (on a single baseline) for a given data 
set. The analyst then may divide the data set into several 
time intervals; residuals in each interval of the experiment 
are then fit with a clock polynomial of their own. Clock 
polynomial analysis becomes very subjective and time- 
consuming. The value of a Kalman f i 1 ter i ng/ smoothing 
approach is thus clear; one does not have to divide a given 
data set into small time intervals for analysis. 

The statistics of the delay residuals are found when 
the clock polynomials are included in the SOLVE least squares 
solution. These statistics are useful in judging the 
goodness of clock (and general) fit of a solution. Least 
squares data analysis continues until the data fit is good; 
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i.e. approaches levels consistent with contemporary data. 

I now present the parameters of particular interest to 
this work, and how the measurement partial derivatives of 
these parameters are calculated for use in the Kalman filter/ 
smoother. The partial derivatives are the components of M(k) 
in equation (III. 4). The parameters of primary interest in 
this research are: the components of earth wobble, x and 
y (mj, m 2 are freely interchanged with x,-y herein); 
the change in length of day, (UTC - UT1 ) ; the Kalman filter 
clock offset and rate terms, x c and x cr ; the atmos- 
pheric wet zenith path delay parameters, L z , and adjust- 
ments to the nutations in obliquity and longitude, A€ and Af 
The clock coefficients and atmospheric parameters are 
clearly "nuisance" parameters in a study of polar motion, 
but the clock and atmospheric delays must be accurately 
modelled to give reliable and precise earth orientation 
results. 

Partial derivatives for polar motion and earth rotation 
are formulated by combining equations (IV. 1) and (A.l) 

(Ma, 1978) 

* b 2000’ s 2000 “(PNSWB terr )»S2000 

Tg = * (IV. 18) 

c c 

where the 2000 and terr subscripts denote the J2000 and 
terrestrial reference frames, respectively. The source unit 
vector S is in the J2000 reference frame. (The spin matrix, 
S, and the source vector, S 20 oo> are not the sarae 
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entity.) The parti a 1 s of the delay with respect to the 
polar motion parameters are found by differentiating 
equation (IV. 18) with respect to x BIH and y BIH ( as 
modified from Lundqvist, 1984) 


9 1 / 9 W \ 

- ( l/c ) ( PNS B terr ] . S 2000 (IV. 19) 

9 ( x BIH>yBIll) \ 9 ( X B 1 H • yB I H ) / 

The partials of the wobble matrix W with respect to x B j H 
and ygjH are (Ma, 1978) 

9W 3Hy(xgiH) 

= R X (yBIH) (IV. 20) 

9x BIH 9x BIH 

and 


3 W 9H x (y BIH ) 

= R y (x BIH ). (IV. 21) 

9yBIH 9 yBIH 

The partial derivatives of the nutation matrix with 
respect to the nutation in obliquity and longitude can be 
readily arrived at using equation (A. 5) 


N = R x (-6)R Z (-A^)R X < 6+Ae ) (IV. 22) 

and differentiating IV. 22 with respect to the obliquity 
angle 


3N 3R X ( +A£ ) 

- R X ( E)R Z (-Af) . (IV. 23) 

3 (AG ) 3(A€) 
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where 


R X (E+AE) = 0 cos(E+AE) sin(£+A£) (IV. 24) 

10 -sin(E+AE) cosjE+AE), 


Then 


9 (AE) 


= R x ( -E ) R 7 ( -A4 1 ) 0 -sin(£+AE) cos(E+AE) (IV. 25) 

,0 - c o s ( £ +AE ) - s i n ( £ + AE ) 


We proceed similarly for the nutation in longitude 


a n aRzi-A^) 

= R x<-fc> fl ~T - H x (£ + A£), 

a(A^) 


(IV. 26 ) 


where 


/ cost-A^) sint-A^) 0 

R z (-A 4 *) = -sin(-A'l') cos ( -A+) 0 


( IV . 27 ) 


“sinA'JJ -cosA^ 0 

= R x (-£) cosAijj “sinA 1 ! 1 0 R x (£+A£). (IV. 28) 


The partial derivative of the delay with respect to the 
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length of day parameter is approximated in a manner similar 
to the wobble partials 


8 x 


3( UTC-UT1 ) 


( 1/c ) 


PN 


as 


8(UTC-UT1 ) 


w ^terr 


) 


S 2000 ( I V . 29 ) 


The partial derivative of S with respect to 
UTC - UT1 , is given by ( Ma , 1978) 

/- sin(GAST) cos(-GAST) 


as 


-cos ( CAST ) -sin(-GAST) 


3(UTC-UT1 ) 




0 


0 


the time of f set , 


0 

0 

0 / 


d(GMST) 
d ( UT ) 
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where GAST is Greenwich Apparent Sidereal Time. 

Atmospheric and clock partial derivatives are found 
differently than the earth orientation measurement partials. 
The effect of the atmospheric path delay, t A , on the total 
delay (t) of a specific baseline is (Ma, 1978) 


L(2) - L ( 1 ) 

t A = (IV. 31) 

c 


where L(l) and L(2) are atmospheric path lengths at VI. B1 
stations (1) and (2). Recall that station (1) is closer in 
distance to the source under observation than station (2) 
at the time the delay is measured. Ma (1978) gives the 
partial of VLHI delay with respect to atmospheric path 
length in his dissertation. 

The partial derivative of delay with respect to the clock 
offset parameter may be difficult to understand at this moment. 
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The meaning of the model should become clear later on. 

The observed delay t includes clock offset effects, tocag. 
Then the measurement partial is simply 

at 

= ±1, (IV. 32) 

8a 0 

where the sign is determined as in the subprogram MARTL (a 
modified form of the subroutine PARTL) within the analysis 
program SOLVE. The expression (IV. 32) is the partial deriv- 
ative for use with the Kalman filter/smoother, and is also 
the partial used in estimating the clock offset constant, 
the first term in the least squares clock polynomial (IV. 17). 
The Kalman filter clock rate partial is assumed to be zero; 
this is an overly simple clock rate partial which models the 
clock ramp to be zero during a VI.BI measurement. Such an 
assumption is only valid (to a fair degree) because VLBI 
measurements on a given source take place over a short 
time (roughly five to ten minutes). The primary partial 
derivatives used herein for Kalman filter VLBI analysis 
have thus been explained. 
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Chapter V. 


EARTHQUAKES AND THE CHANDLER WOBBLE (A REVIEW) 

The main thrust of this work is to clarify whether 
earthquakes are a significant excitation source of the 
Chandler wobble of the earth. Apparently, other mechanisms 
may be able to drive the wobble, such as atmospheric mass 
movements or core mantle coupling, but the literature rele- 
vant to these phenomena will not be widely explored here 
(see Wilson et a 1 . , 1976 and Barnes et a 1 . , 1983). Liter - 
ature concerning earthquake excitation of the Chandler 
wobble will be examined. 

Interest in the effects o f earthquakes on polar motion 
in the late 1960’s led to the organization of a NATO 
Advanced Study Institute ori "Earthquake Displacement Fields 
and the Rotation of the Earth" (Mans inha, S m y 1 i e and Beck, 
1970 ) , which was held in 1969 . The edited volume based on 
this conference provides a good historical basis for any 
earthquake excitation polar motion study. Prior to 1965, 
the changes in the earth’s inertia tensor due to earthquakes 
were thought to be inadequate to excite Chandler wobble. 
Press (1965) demonstrated that mass displacements during a 
great earthquake extend over distances much further than 
were thought previously. Perhaps, since more mass is 
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involved in earthquakes than was previously thought, 
wobble may be excited by the larger associated inertia 
tensor changes. 

An example of the involvement of mass at distance from 

a large earthquake and its effect on Chandler wobble 

is given by Stacey (1977). He considers an exemplary 

inertia tensor change for maintaining the Chandler wobble to 

2 8 2 

be about 1.3 x 10 kg-m , and considers that great 
earthquakes, such as that of 1964 in Alaska, may be the prime 
movers of the Chandler wobble. The inertia tensor change 
due to an Alaskan-like earthquake is found (Stacey, 1977) 
by approximating the earthquake geometry by two blocks of 

3 

dimensions 800 km x 200 km x 200 km, and density 3000 kg/m , 

separated by a mean distance of 200 km, and displaced by 

22 meters. The resulting inertia tensor change computed 

2 8 2 

using the simple earthquake model is 8.5 x 10 kg-m , 

which obviously falls short of Stacey's excitation require- 

2 8 2 

ment of 1.3 x 10 kg-m . 

However, if mass displacements at larger distances 

are included in the calculation of the perturbation of the 

inertia, Stacey (1977) states that the inertia tensor 

change due to the earthquake and associated mass movements 

is on the order of the required moment adjustment, 1.3 x 
2 8 2 

10 kg~m . Thereby a great earthquake event might 
be able to significantly excite the Chandler wobble. Still, 
many geophysicists are not convinced that a great earthquake, 
occuring once or twice per decade, is sufficient to maintain 
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the Chandler wobble, though even they are not very sure of 
the time required for the Chandler wobble to damp out 
(Stacey, 1977). 

More compelling evidence in support of the forcing of 
Chandler wobble by earthquakes was given by Mansinha and 
SmyJie (1970) at the NATO Advanced Study Institute. Recall 
that Figure 1 depicts the modification of the Chandler 
wobble pole path due to a sudden mass shift in the earth. 

The rotation pole orbits about a shifted secular pole after 
excitation, with the orbital radius generally of different 
magnitude than prior to excitation. Mansinha and Smylie 
(1970) exploit the fact that a "kink" (change in curvature) 
is generated in the polar motion path due to a mass shift and 
try to correlate kinks in real data with large; earthquakes, 
of Magnitude M > 7.5. The polar motion data (with annual 
wobble component removed) were from the BIH and ILS-1PMS 
(International Polar Motion Service) for the years 1957.0 to 
1968 . 0 . 

The procedure by which kinks in the data are revealed 
is now described. If polar motion is continuous, the 
rotation pole moves through angle a in 6 days with 
respect to the "mean pole", with angle a given by (Mansinha 
and Smyl ie , 1970 ) 

2k 

a = x 6 , ( V . 1 ) 

438 

where 438 is their approximate period i ri days of the Chandler 
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wobble. A kink due to a shift in mass will disrupt the 
observed continuous polar motion. Mansinha and Smylie 
(1970) start their analysis, say at 1957.0, and fit a 
curve between two pole positions of the BIH data set (the 
analysis is done independently for ILS-IPMS data). ILS-1PMS 
and BIH measurements are made at regular time intervals. 

After the first arc is constructed between the two points, 
the arc is extended by assuming the same circular path 
(a predictive step) to approximate the next BIH datum. The 
predicted datum location on an x- y polar motion plot is 
then compared with the actual BIH datum location. If the 
prediction falls within some acceptance distance, a, of the 
actual BIH datum, no "kink" in the data has been observed, 
and a new arc is fitted to the three BIH data points in order 
to continue the comparison procedure between actual and 
predicted data (Mansinha and Smylie, 1970). The acceptance 
distance, a, is proportional to the error in the BIH measure- 
ments (a has centi-arc second units). 

If the actual BIH polar motion position is further from 
the predicted datum than the acceptance distance "a", a break 
or "kink” in the pole path has been detected (Mansinha and 
Smylie, 1970). A new arc is now created using the aberrant 
BIH datum as a starting point and the comparison process 
continues again. The breaks found in the ILS-IPMS and BIH 
data were then compared against (M > 7.5) earthquake event 
times for 1957.0 to 1968.0 (Mansinha and Smylie, 1970), and 
random probability criteria were used to determine if the 
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breaks and the earthquakes were correlated. "For a = 0.015, 

12 out of 22 earthquakes fall within tlO days of a break and 

17 out of 22 earthquakes fail within ±20 days of a break. 

The corresponding values of RP [random probability) are 
~3 - 3 

1.3 x 10 and 1.6 x 10 respectively." (Mansinha and 
Smylie, 1970). The random probability function used in the 
study (cumulative binomial distribution) was (Mansinha and 
Smy) ie , 1970 ) 


RP 


n / n 

A l; 




( V . 2 ) 


"The elementary probability p is the probability of a date, 
drawn at random, falling within co [to = correlation window; 

= 26 ] of a break [kinkj and is given by the proportion 
of time axis of the total span of 11 years that is occupied 
by segments w wide on each side of a break." (Mansinha and 
Smylie, 1970). Results are for BIH data sampled every 
ten days. Considering that only 22 large events occurred 
during the 11 year data span, and that 17 out of 22 
earthquakes fall within ±20 days of a break, one might be 
inclined to believe that the polar motion breaks and 
earthquakes are well correlated. 

Unfortunately, the study by Mansinha and Smylie (1970) 
may have been ahead of its time, although the method has a 
simple and believable mathematical basis. The polar motion 
data used in the analysis was only sampled at ten-day inter- 
vals, which leaves much to be desired in time resolution. 
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Also, Mansinha and Smylie (1970) note that the uncertainty 
of the polar motion data was on the order of ten milliarc- 
seconds , which is much larger than present-day uncertainties. 
However, current VLBI polar motion data sets do not 
approach the length (11 years) of the BIH and ILS-IPMS data 
used by Mansinha and Smylie. As a final point, Morabito 
and Eubanks (1985) intercompare ILS-IPMS polar motion data 
with concurrent polar motion component estimates determined 
by VI. B I , satellite laser ranging and lunar laser ranging 
techniques and find differences between modern data and ILS- 
IPMS data of 50.3 mi 1 1 iarcseconds rms in the X component of 
wobbble and 23.3 m i 1 1 i arcseconds rms in the Y component; 
these differences are primarily due to seasonal variations 
in the ILS-IPMS data, but considerable differences still 
exist even after the seasonal variations are removed from 
the data. Morabito and Eubanks (1985) state that "These 
results indicate that any geophysical conclusions derived 
from ILS data should be interpreted with caution." The 
intercompared data spanned 1976-1979. Morabito and Eubanks' 
warning should be strongly heeded. Nevertheless, the work 
of Mansinha and Smylie (1970) was preeminent, and acted as 
a first step in contemporary Chandler wobble excitation 
research. 

Dahlen (1973) uses a SNKEI ("... spherically symmetric 
non -rotating Earth which is in hydrostatic equilibrium and 
has an isotropic elastic constitutive relation."- Smith, 

1974) Earth model to calculate the pole shifts due to the 
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1960 Chilean and 1964 Alaskan earthquakes. 


T h o changes In 


the product of inertia components C23. C 23 (which are of 

primary interest in polar motion excitation see equations 

11.7 and 11.9) of the earth model due to an earthquake at 
colatitude 0 O and longitude 4> 0 are (Dahlen, 1973) 


AC j 3 =Mq{ 1'2 (h) [ [sin2asin6cos\t-l/2cos2asin26sinX.) 
xsin20ocos<Do'-2(l/2sin2asin26sin\-cos2«sin6cos\) (V.3) 

xsin0osin®o] +r 2^)[ _s;ln2 ^ s ^ n ^ s ^ n2 ®O cos<I) o] 

+ F 3 ( h ) [ (sinaoos26sinX.-cosacos6cosX.)eos20QCOs<J>Q 

+ ( sinacos6cos\+cosacos2{is.inX.) cos0()S.in<I>o] ) 

AC 2 3 = Mg { T 2 (h) [ (sin2osin6cos\<-l/2cos2asin26sin\) 
x sin20os.in<I»o + 2fl/2sin2a8in26sinX-cos2asin8cos\) (V.4) 

.'■v 

x s 1 n0ocos®o + r 2 ( h ) [-sin26sin\sin20osin<l>o] 

+ 1’ 3 ( h ) [sinaeos26sin \ - cosacos6cos\)cos20QS.in<J>o 
-(sinacos6cos\. + cosaeos268inX.)eo80gco8<l>g|, 

where h is the earthquake depth, 6 is earthquake; dip, a is 
the strike, \is the slip angle and Mg is the scalar moment 
of the earthquake. r 2 ( h ) , [’2(h), and r 3 ( h ) are canonical 
functions of depth computed for the SNREI Earth model under 
consideration; details about them can be found in Dahlen 
(1973). Clearly the product of inertia perturbations AC23 
and AC 2 3 of equations (V.3) and (V.4) are directly 
proportional to the seism. it: moment Mg of a given 
earthquake; thus AC 2 3 and AC 23 are critically 
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dependent on moment determinations. 

The pole shifts due to ACjg and AC 23 are calculated by 
Dahlen (1973) for the i960 Chilean and 1964 Alaskan earth- 
quakes, and are compared with pole path shifts as seen in 
BIH polar motion data (Smy.lie and Mansinha, 1968). Dahlen 
(1973) finds that his theoretical polar motion shifts are an 
order of magnitude less than those observed by Smylie and 
Mansinha (1968), although the pole path shift directions 
determined with each method are in agreement. Noise contam- 
ination is blamed by Dahlen (1973) as the phenomenon which 
precludes observation of large polar motion path shifts in 
BIH data. He also concludes that "The direct observation 
of the effect of a large earthquake on the path of the 
Earth's rotation pole will not be possible until some way 
is found to reduce the noise contamination in the 
observed polar motion path by at least a factor of ten." 

VLBI data is more precise than the BIH data of 1968 by at 
least a factor of five to ten, so kinks in polar motion 
paths due to large earthquakes might be observed presently. 

Dahlen (1973) next examines the cumulative effect of 
many large earthquakes on Chandler wobble power P, which 
is computed via the formula 

Q 2 

P = <R > ( V . 5 ) 

u> o 

where Q is the Chandler wobble quality factor, cjq is the 
characteristic oscillation frequency of the wobble, and 


73 


2 

<R > "... is the mean squared polar shift, per unit time 

2 

associated with the earthquake sequence.” Since Q and <R > 
are not accurately known, any power estimated using (V.5) 
is only of limited accuracy. D a h 1 e n accounts for this in 
some fashion by considering several values of Q in his 
power calculations, and he also compares the magnitude 
moment relations of Brune (1968), and Aki (1967, 1972) 

and the effect of each on the power calculations. The 
magnitude-moment, relations are displayed in Figure 7 (from 
Dahlen, 1973). 

2 

It appears that the Aki u> model (1967, 1972) provides 

a better fit to observed surface wave and field data than 

Brune' s (1968) magnitude -moment relationship, if the Sanriku 

(1933) datum is not included in the fitting process. There 

is some justification for not including the Sanriku 

earthquake in forming a general magnitude moment relation. 

O'Connell and Dziewonski (1976) state that it is difficult 

to assign a source mechanism to Sanriku since the event 

occurred near the junction of the Eurasian, Pacific and 

Philippine plates, and its moment is somewhat low for large, 

shallow, subduct: ion-zone events. In any event, use of Aki’s 
2 

(1967, 19 72) co moment -magnitude relation gives a more 

favorable outlook (by power comparisons) of earthquake exci- 
tation of the Chandler wobble than does the Brune (1968) 
relationship . 

For the calculated Chandler p o w e r values as a function 
of Q and moment estimation process, please refer to Dahlen 
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Fig. 7. Seismic moments M 0 vs. earthquake 
magnitude, as proposed by Brune (1968) and Aki 
(1967, 1972). (From Aki, 1972 and Dahlen, 1973). 
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(19 7 3). IJ a h 1 e n (1973) ultimately d c c i d e s "... that the 
Earth's seismic activity can account for no more than 10 per- 
cent (for Q - 100), and probably even significantly less, of 

the observed level of excitation of the Earth's Chandler 
wobble." But since the value of the quality factor, Q, is 
ill-constrained, it is difficult for anyone to know if such a 
conclusion is valid, or certainly final. 

In a much different vein, evidence was found by Kanamori 
and Cipar (1974) to support the hypothesis that earthquake 
events are accompanied by large, slow, aseismic deformations. 
The foundation for their study comes from a long-period 
strain seismogram recorded during the Chilean earthquake of 
22 May, 1960. On the seismogram, accompanying the P wave 
arrival from a M s = 6.8 foreshock, which occurred fifteen 
minutes prior to the main Chilean event, Kanamori and Cipar 
(1974) observe a long-period wave of period 300-600 seconds. 
That such a signal is real may be a valid question, but 
Kanamori and Cipar (1974) respond to this by stating that 
the seismogram shows a quiet, noise free trace prior to the 
fores hock, and also note that such a long period precursor 
may be evident for the Chilean 1960 event, but not for other 
earthquakes because the Chilean earthquake is one of the 
largest events ever recorded using "modern” seismographs. 

They do admit, however, that the long-period precursor signal 
could be due to instrumental instability. 

Next, Kanamori and Cipar (1974) determine fault param- 
eters for the 1960 main Chilean earthquake and the long- 
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period precursor by matching synthetic seismograms to the 

observed strain seismogram traces. They find a moment of 

30 30 

2.7 x 10 dyne-cm for the main shock and 3.5 x 10 dyne-cm 

for the aseismic long-period precursor. The aseistnic moment 

is of the same order as the moment of the elastic main shock. 

Kanamori and Cipar (1974) postulate that the Chilean sequence 

occurred as follows: "About fifteen minutes before the main 

shock a gradual slip having time constant of about 450 sec 

started taking place at the boundary between the downgoing 

oceanic lithosphere and the relatively weak as thenosphere 

beneath the [South American] continent. This slip caused a 

stress concentration at the lithosphere lithosphere boundary 

above. The large foreshock may represent a brittle fracture 

caused by such a stress concentration. Eventually this 

stress concentration exceeded the breaking strength of the 

1 i t hos phe r e - 1 i t hos phe r e boundary causing a major catastropic 

failure, the main shock." They also suggest that the Chilean 

aseismic slip may have occurred for hours or days before the 

main shock and feel that precursory long-period phenomena, 

even if only observed ten to fifteen minutes before a main 

shock, may be useful in preparing human populations near the 

epicenter for an impending event. Obviously, further study 

of aseismic stress releases is warranted. 

Press and Briggs (1975) find correlations between 
Chandler wobble and earthquake activity through a pattern 
recognition analysis, the details of which are too many to 
be included here. They conclude that earthquakes of certain 
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positions (latitude, longitude, slip and fault orientation) 
are associated with years of increasing Chandler wobble 
amplitude, and hypothesize that large preseismic and post- 
seismic deformations (over hours or days) play a part in 
Chandler wobble forcing, in addition to the slip occurring 
during the main earthquake shocks. Press and Briggs (1975) 
also note that earthquakes along the Kurile -Japan- Marianas 
seismic belt and in China do not show a correlation with 
increasing Chandler wobble amplitude, and suppose that 
" . . . ; perhaps slow deformation occurs in these belts during 
aseismic periods many years before or after major 
earthquakes . " For the most part, Press and Briggs (1975) 
conclude that Chandler wobble excitation is due to large 
mass displacements associated with earthquakes. 

Synthetic curves which mimic ILS polar motion data (with 
annual wobble removed) have been generated by O'Connell 
and Dziewonski (1976) for the time period 1901-1970, by the 
superposition of polar motion shifts from 234 large* earthquakes 
of M s >7.8. The polar motion shifts are found from 
changes in the earth inertia tensor, ACjj, which are related 
to the earthquake moment tensor via the following equa- 

tion (O'Connell and Dziewonski, 1976) 

AC i j =0ijkl M kl* (V.6) 

where is a tensor with components determined by 

choice of earth model, and latitude and longitude of a 
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given earthquake epicenter. For details, see O'Connell 
and Dziewonski (1976). The source moments M k i for the 234 
earthquakes are approximated using simple double couple 
moments (of scalar moment Mq) in unison with plate 
tectonic information (O'Connell and Dziewonski, 1976). 

The seismic moments Mq are calculated from surface 
wave magnitudes M s reported by Duda (1965), and other 
sources, using the relation (O'Connell and Dziewonski, 

1976 ) 


Log 10 M 0 = 8 8 + 2 5 M s • (V.7) 

Once the pole shifts are calculated, the synthetic wobble 
m(t) is found with the expression (O'Connell and Dziewonski, 
1976 ) 

m ( t ) =m ( t 0 )exp{ i u( t -t 0 ) ) + (V.8) 

+ 2 1 s k H(t-t k ) {1- (l+w 0 /Q)exp|iw(t- t k ) ] } , 
k 

where m(tg) is the reference (starting) position of the 
Chandler wobble pole at time equals tg # Cx) ^ 0 (1 - i/2Q) . 

wg is the Chandler wobble frequency, Q is the Chandler wobble 
quality factor, H ( t - 1 k ) is a step function, s k is the pole 
shift due to an earthquake at the epoch t k , 0 is the diurnal 
rotation frequency, and t is the time for which the synthetic 
signal m(t) is produced. The real and imaginary parts of 
m(t) are the Chandler wobble components, Xj and -X 2 respec - 
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t i ve 1 y . 

Synthetic polar motion curves for the Xj wobble component 
(for various Q values) and smoothed ILS data with annual 
wobble removed are presented in Figure 8. The curve 
for Q s =-• Qsynthetic r= 100 resembles the ILS curve in general 
form* but the match is by no means exact. O’Connell and 
Dziewonski (1976) then compute power spectra of the various 
synthetic and ILS curves and compare the Chandler wobble 
peaks. In general, the agreement of ILS and synthetic data 
is fair. O'Connell and Dziewonski (1976) conclude that 
earthquake-associated mass movements may account for about 
fifty per cent, of Chandler wobble excitation . 

There is a major problem with the work of O'Connell 
and Dziewonski (1976). Kanamori (1976) points out that the 
surface wave magnitude used in Duda (1965) and by O'Connell 
and Dziewonski (1976) is not the same as the M s referred to 
in the relation Logjo Wq = 8.8 + 2.5 M s . The net effect of 
this error is that O'Connell and Dziewonski (1976) overesti- 
mate the scalar moment Mg of many of the earthquakes used in 
their study by a factor of 5.6 (Kanamori, 1976), and thereby 
overestimate many of the Chandler wobble shifts due to these 
earthquakes. This casts doubt upon the results of O'Connell 
and Dziewonski (1976). Another problem with O'Connell and 
Dziewonski 's work is that their comparison was made in the 
time domain, as opposed to the excitation domain; Chao 
(1985) addresses this problem in detail. Kanamori (1976) 
also states that the direct effect of earthquakes on Chandler 
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Fig. 8. Simulated polar motion curves resulting 
from excitation by large earthquakes; the damping 
for each is determined by Q s . ILS is smoothed 
observed polar motion with annual term removed. 
Only the component along the Greenwich meridian 
is shown. (From O'Connell and Dziewonski , 1976). 
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wobble as calculated In a correct manner (using proper M s 
values) is probably an order of magnitude too small to give 
observed wobble amplitudes. He does acknowledge that large 
aseismic deformations accompanying earthquakes may indeed 
be responsible for driving the Chandler wobble, and proposes 
that these deformations may be detected in the future with 
global arrays of ultra long period se i smo 1 ogi ca 1 instruments. 

In his 1977 article, Kanamori asserts that the Gutenberg 
and Richter (Gutenberg, 1956) magnitude energy relation 
log E = 1 . 5M + 11.8 is of questionable applicability to 
earthquakes of rupture length of 100 km or more (but works 
well for smaller rupture lengths). E is the energy released 
from an earthquake and M is the magnitude determined from 
seismic records. "This [problem! arises from the fact that 
for such a great earthquake the magnitude M which is deter- 
mined at the period of 20 s (or converted from m (body wave 
magnitude) determined at shorter periods) does not represent 
the entire rupture process of an earthquake. In fact, 
there is little correlation between M and the rupture length 
for great earthquakes." (Kanamori, 1977). 

Kanamori then derives, following laborious and rigorous 
analysis of large earthquakes, an alternate energy -magnitude 
relation for great events 

log W 0 = 1 .5 M w + 11.8. ( V . 9 ) 

M w is the new magnitude for great earthquakes and W 0 
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is the minimum strain energy drop for great earthquakes. 
According to Kanamori (1977), the minimum strain energy 
drop is equal to the seismic wave energy. Also he states 
that W 0 can be calculated using the approximation 

w 0 * M 0 / ( 2 xl 0 4 ) . (V . 10 ) 

where Mg is the seismic moment, of an earthquake (CGS 
units) . 

The minimum strain energy drop W 0 (and seismic wave 
energy released) is then calculated for (shallow) great 
earthquakes from 1904 to 1976 (Kanamori. 1977) and the 
results (5-year running averages) as a function of time are 
presented in Figure 9. In addition to Wq. the Chandler 
wobble amplitude; N, the annual number of earthquakes of 
M s > 7.0 (a five-year running average); and a five- 
year running average of E (explained earlier) are included 
in the figure. The correlation between Wq and wobble 
amplitude is high, and it appears that the wobble amplitude 
leads Wq. N seems to lead the Chandler wobble ampli- 
tude. 

Kanamori (1977) proposes one explanation for the 
results in Figure 9. "One possible mechanism that accounts 
for the correlation between the wobble, Wq. and the 
activity of moderate to large earthquakes (as indicated by 
N | is that an increase in wobble amplitude triggers 
worldwide seismic activity and accelerates plate motion, 
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Fig. 9. Correlation between the amplitude 
(envelope) of the Chandler wobble, W 0 (5 year 
running average), annual number N of earthquakes 
of M s > 7.0 (5-year running average), and 
the Gutenberg Richter energy E (5-year running 
average). (from Kanamori, 1977). 





which eventually leads to great decoupling earthquakes [i.e. 
decoupling along major plate boundaries]. This decoupling 
causes the decline of moderate to large earthquake 
activity. " 

Kanamori's (1977) analysis is straightforward. An 
alternate magnitude scale to describe the magnitude-energy 
relationship of great earthquakes is necessary. The 
correlation of wobble amplitude and energy released by great 
earthquakes is striking. The data do indicate that a 
Chandler wobble amplitude peak does precede a maximum in 
the seismic wave energy (1950 to 1965). The paper does 
not include the effects of aseismic slip. Kanamori's 
(1977) work suggests that the Chandler wobble triggers 
earthquakes, as opposed to vice versa. 

It was mentioned earlier in this work that the Chandler 
wobble is a long period normal mode of the earth. It is 
because of this long periodicity that the Chandler wobble is 
of geophysical interest; few geophysical phenomena have been 
observed (or exist) at this frequency. Smith (1977) looks 
at the excitation of the Chandler wobble by seismic events 
by treating wobble as a normal mode. He uses his (1974) 
equations describing "... the infinitesimal free elastic- 
gravitational oscillations of a rotating, slightly elliptical 
Earth ... " to compute polar shifts due to the I960 Chilean 

(modelled as a two-event source) and 1964 Alaskan earthquakes, 
which are two of the largest seismic events in recent history. 
The details of Smith (1974) are not fully described here. 
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Smith's (1977) normal mode calculation of Chandler wobble 
shifts differs from elastic dislocation methods typically 
applied to the excitation problem; it also avoids difficul- 
ties encountered in equations of motion describing the static 
deformation of the fluid core of an earth which is not 
rotating, which are used in such elastic dislocation calcu- 
lations (Smith, 1977). Thus, Smith's (1977) analysis might 
be used to determine which of the many quasi-static Chandler 
wobble excitation studies is most correct. 

Smith (1977) represents the particle displacement field, 
s(r), of the Chandler wobble normal mode by a truncated 
series of spheroidal vector fields, oj™ , and toroidal vector 
fields, t i™, 

s ( r ) = + o 2 ±1 (r) + t 3 ±:l (r) (V.U) 

where r is a radius variable. Equation (V.ll) is inserted 
into the differential equations of motion of Smith (1974), 
and the solutions of these equations approximate normal mode 
eigenfunctions describing the Chandler wobble. The excita- 
tion theory developed by Dahlen and Smith (1975) is then 
applied to the calculated eigenfunctions to find the 
Chandler wobble shifts due to the aforementioned 
earthquakes (Smith, 1977). 

Smith (1977) argues that the truncated series (equation 
V.ll) is accurate enough for his Chandler wobble study. He 
bases this on the fact that "... the principle response of 
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an elastic solid body to rotational disturbances was charac- 


terized by 1=2 spheroidal distortion." (Smith, 1977). Also 
Smith (1977) suggests use of equation ( V . 1 1 ) is valid 
because (to first order in ellipticity) equation (V.1J) is 
of the correct, form "... which represents the response of 
an incompressible, homogeneous, inviscid rotating fluid core 
to slight rotations of a rigid mantle...." Equation ( V . 1 1 ) 
does include an 1=2 spheroidal term. However, it is 
p r o b a b .1 y m o r e r e a .1 i s t i c to c o ri c .1 u d c that Smith used the 
truncated particle displacement series < V . 1 1 ) so that the 
computations and analysis could be carried out with reason- 
able effort. Smith's (1977) calculations should be 
performed again using more terms in the particle displace- 
ment series. 

Smith (1977) presents the polar motion shift results 
in terms of a magnitude and a direction, in degrees East of 
Greenwich, rather than by two wobble components. For 
comparison, values found by Dahlen (1973) using a quasi 
static calculation are given in parentheses. For the first 
event of the 1960 Chile source model (Chile (1)) the 
magnitude is 0 0 2 1 1 6 (070210) and the d i rection is 114.3° 
(116°). For the second event (Chile (2)) the values are 
0702803 (070280) and 117.7° (118°). Finally, the shift 
magnitude and direction for the 1964 Alaska earthquake are 
0700723 (070073) and 201.3° ( 202°) : (Smith, 1977 ). 

Obviously the agreement of the polar shifts and 
associated directions between Smith (1977) and Dahlen (1973) 
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is very good, implying that the quasi-static calculation of 


Dahlen (1973) is among the better attempts to determine polar 
shift magnitudes using a dislocation source on a non-rotating 
earth. No conclusion is made in Smith's (1977) paper as to 
whether earthquakes are the primary source of Chandler 
wobble excitation. 

Marisinha, Smylie and Chapman (1979) revisit the Chandler 
wobble, and again conclude that earthquakes are the major 
driving mechanism of the Chandler wobble. They correct the 
static dislocation theory of Smylie and Mansinha (1971) for 
errors made in transforming planar dislocation models to a 
more realistic spherical geometry. These corrections bring 
pole path shift amplitudes and directions (degrees East of 
Greenwich) of Smylie and Mansinha (1971) into better 
agreement with the estimates of Smith (1977), with the 
resulting amplitudes being different by no more than a 
factor of one-half to two. Mansinha, Smylie and Chapman 
(1979) use the SNREI model of Dahlen (1971, 1973) in their 

computations of the shifts, and speculate that pole shifts 
observed (astronomically) in the future may be used to 
constrain focal mechanisms of earthquakes. 

The cumulative effects of earthquakes on wobble is then 
studied by Mansinha, Smylie and Chapman (1979). They 
characterize the root -mean-square (rms) Chandler wobble 
amplitude using the relation 
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{ “ ( s o / s m ) 
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b/ c 


N 0 t 


2 


} 1/2 , ( V . 12) 


" . . . where s m is the average polar shift due to the largest 
magnitude considered, sq is that due to the smallest magni- 
tude considered (with N q the c o rrespo n ding f r equoncy ) , and 
t is the damping time of the Chandler wobble.” "b" and ” c " 
are constants derived from frequency -magnitude and pole path 
shift -magnitude relations respectively. Mansinha, Smylie 
and Chapman (1979) then use the 30 largest earthquakes from 
1901 to 1964 as taken from O’Connell and Dziewonski (1976), to 
form a plot (Figure 10) of rms Chandler wobble excitation 
as a function of the constant " c ” mentioned previously. They 
assume a damping time of 20 years (Q ; 52). 

Mansinha, Smylie; and Chapman's ( 1979) calculations as 
shown on Figure 10 surely support the view that earthquakes 
are, at worst, a large contributor to C h a n d I e r wobble 
excitation and, at best, the primary excitation source. 
Mansinha, Smylie and Chapman (1979) continue to use BIH data 
as the observational standard of comparison for their 
theoretical results, and base; cumulative earthquake excita- 
tion conclusions on their magnitude-moment of choice. 

Research which indirectly supports the hypothesis that 
earthquakes primarily excite the Chandler wobble was done by 
Wahr (1983) in an attempt to determine the level at which 
atmospheric phenomena drive the Chandler wobble. Wahr 
(1983) considers earthquakes and the atmosphere as the most 
plausible means of exciting the Chandler wobble, but 
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Fig. 10. Cumulative Chandler wobble excitation by 
earthquakes. Shaded area indicates domain of 
possible excitations (boundaries are extremes). 
Range is from 0.08 to 0.50 while observed level 
is 0.15. Some of the more extreme values of 
seismic moment would provide too much Chandler 
excitation. S ^ is a reference pole shift 
associated with a reference magnitude Mj; 
see reference for details. (From Mansinha, 

Smylie and Chapman, 1979). 


concludes that the atmosphere is not the main source of 
excitation. Atmospheric data and models are used by Wahr 
to compute power, coherence and phase spectra of atmospheric 
excitation sources such as atmospheric pressure, the 
pressure- driven ocean, mountain torques due to variations in 
pressure at different locations on a mountain, and wind drag 
over the surface of the earth. These spectra are compared 
to excitation spectra from ILS data for the years 1900-1973 
to investigate if atmospheric effects do drive Chandler 
wobble. 

Wahr (1983) considers several ways to look at mechanisms 
driving Chandler wobble: one way is by examining excitation 
functions, ^ . and associated derived quantities such as 
excitation power spectra in frequency space; another is by 
examining the wobble components (m i( m 2 ) in the time domain. 
Wahr uses the former technique. The ultimate result of 
Wahr's study is that no more than 20-25 percent of the 
Chandler wobble observed (ILS) excitation is due to atmos- 
pheric and oceanic sources. Wahr (1983) does not take into 
account the effects of ground water storage or El Nino 
events . 

In more recent work. Gross and Chao (1984) deconvolve 
LAGEOS (satellite for laser ranging) Chandler wobble polar 
motion data for the period 2.1 January, 1977 to 20 January, 
1984 to look at the excitation function of the wobble. They 
observe a shift in the excitation axis by 0'.’040 in the -54 
degrees East longitude direction during the time period 


8/13/77 to 9/2/77, and assume the shift is due to the Sumba 


earthquake of 8/19/77. The correlation in time of the two 
events is quite good. However, when Gross and Chao (1984) 
compare their observation with Dahlen's theoretically 
derived value of the Sumba shift due to elastic dislocation 
models (O'.'OOOS towards 93 degrees East longitude), they 
they find little agreement between theory and experiment. 
Gross and Chao ( 1984) conclude that the 1.AGE0S observation 
may have resulted from some discontinuity in processing 
LAGEOS data. 

Gross and Chao (1984,1985) also consider the possibility 
that the descending slab associated with the 1977 Sumba 
earthquake may have broken away from its parent oceanic 
lithosphere and excited the Chandler wobble in its downgoing 
motion; but they reach no conclusion about this hypothesis. 
They (1985) also look for a correlation between polar motion 
excitation and the 1982-83 El Nino and believe they observe 
one. The 1982-83 El Nino was particularly strong and may 
have been involved in Chandler wobble excitation. 

Gross (1986) compares Chandler wobble excitation 
function power spectra derived from theoretical static 
deformation fields of 1287 large earthquakes during 1977- 
1983 (a time involving one great earthquake) with power 
spectra of observed excitations derived from LAGEOS earth 
orientation measurements. He finds that the earthquake 
(amplitude) power is about 56 deciBels (dB) less than the 
power necessary to drive the observed Chandler wobble. 
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Gross concludes that the static deformation fields of 
earthquakes during 1977-83 are not the primary excitation 
source of the Chandler wobble, but also notes that no ex- 
tremely large earthquakes, such as the 1960 Chilean or 1964 
Alaskan events occurred during this time interval. He 
leaves open the possibility of Chandler wobble excitation 
by aseismic slip, whose contributions to Chandler wobble 
have not been currently evaluated. 

In another investigation, Robertson (1985) examines IRIS 
VI,BI polar motion components plotted against one another 
over five-day intervals during 1984 and early 1985. He 
observes variations in the data similar to the breaks 
observed by Hansinha and Smylie (1970), but is hesitant to 
correlate the variations with another large Chilean 
earthquake which occurred in eariy 1985. Perhaps it is 
better to be cautious in interpreting early VLB1 results. 

Hinnov (1985) takes a much different route to explain 
the Chandler wobble excitation by water storage and air mass 
loading. She uses Northern Hemisphere precipitation and 
temperature observations for the years 1900-1979 to form an 
excitation series. The air mass and water storage excitation 
data display almost zero phase lag and significant coherence 
with ILS polar motion observations near the Chandler 
frequency. Hinnov (1985) believes that the aforementioned 
meteorological phenomena can fully explain Chandler wobble 
excitation. 

It is of interest to note that geophysical observers do 
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do not always see eye-to-eye with theoreticians regarding 
numerical agreement of strain magnitudes. The following are 
two early quotes on the subject. "Surface strain observa- 
tions are at least an order of magnitude greater than 
predicted by the elastic theory of dislocations." (Mansinha 
and Smylie, 1970). And from a question at the 1969 NATO 
Advanced Study Institute about the effects of earthquakes 
on earth rotation: " Although you [Richard Haubrich, Uni- 

versity of California/San Diego] mentioned that the exci- 
tation from historical [earthquake] events is at least an 
order of magnitude too small to maintain the observed wobble, 
1 [M.W. Major, Colorado School of Mines] am not clear as to 

how one arrives at the magnitude of the excitation. The 
main problem is that most of us see strain steps at distances 
which are embarrassingly large compared to those predicted 
by either the half-space or spherical dislocation theory." 
(Mansinha, Smylie and Beck, 1970). 
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Chapter V I . 


KALMAN FILTERING OF VI.BI DELAY DATA FOR EARTH ROTATION AND 
ORIENTATION PARAMETERS 

A Kalman filter has been designed to make optimal 
estimates of earth rotation and orientation parameters 
which are derived from Very Long Baseline Interferometry 
delay data. Generally the filter will be applied to one 
to two day long VLBI data sets. Such a data set will have 
several hundred or more delay observations, and should allow 
one to estimate polar motion and earth rotation variations 
over times of a few hours to two days. If one desires 
observational data of longer duration, one can examine mul- 
tiple data sets . 

Earlier it was mentioned that the derivation of equation 
(III. 5) 

x k + 1 » <M k + 1 , k ) x k + I’ (k ) w k 

would follow, and this will be handled now. Equation (III. 5) 
describes the dynamical behavior of the state vector x k in 
time. Recall <D(k + l,k) is the state transition matrix 
which models deterministic changes in x as one proceeds from 
time t k to time t k + 1 . w k is a stochastic noise input in 
the state equation, with T(k) accounting for proportion- 
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aiity. Equation (111.5), or relations iike it, which 
describe the physical variations of the estimated parameter 
vector x, are usually found from continuous (as opposed to 
discrete) differential equations characterizing the physical 
system . 

The continuous differential equation form of (III. 5) is 
typically given as (in notation of Gelb, 1974) 


x(t) = F ( t ) x ( t ) + G(t)w(t) , 


(VI. 1) 


where x(t) is the state vector, w(t) is the stochastic 
forcing function and F(t), G(t) are proportionality matri- 
ces. Consider, as a first step towards understanding the 
system represented by (Vl.l), that there is no random 
forcing w(t) and that F(t) is time invariant. The differ- 
ential equation representing this simplified system is just 

dx ( t ) 

= F x ( t ) . (VI. 2) 

dt 


If one rearranges the above differential equation and inte- 
grates from time t k to t ktl 


x k + l 

f dx ( t ) 


J x( t ) 

x k 


t k+ 1 

F dt , 
t k 


or one can readily find that 


(VI . 3) 
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x k + l 


l c* x 1 » { F ( t k + 1 t k )}| x k . 


(VI .4 ) 


where x k , x k + ^ are the parameter values at times t k , t k ,i 
respectively. See* also Geib et al. (1974). 

Thus the state vector x k * j is related to the state 
vector at the previous time x k by the quantity exp{F(t k+ j~t k )}. 
The matrix exponential is defined by (Geib et al.. 1974) 

A A 

e A = I -* A + + t (VI. 5) 

2 ! 3 ! 

where I is the identity matrix. The state transition matrix 
® ( k + 1 , k ) is 


®(k+l,k) = exp(F(t k+1 t k )) 


(VI . 6 ) 


T h e state transition matrix depends only on the time interval 
(t k . t -| - t k ) as should be the case for a stationary system. 

Equation ( V I . 6 ) relates the continuous matrix V and the 
discrete state transition matrix, <D ( k M , k ) . 

Through a derivation procedure similar to the above, one 
can find the relationship between state vectors x k , x k+ j for 
a non -stationary process with no random forcing 


x kH 


[ { e x p ( 


i k + l 

F(t) d t ) > ] x k 
^ k 


(VI . 7 ) 


Evidently, the state transition matrix now takes the form 
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0(k+l ,k) - exp ( 


t k+ 1 

F(t) dt) . 
t k 


(VI .8) 


Using the general form (VI. 8) for $ ( k -t 1 , k ) one can readily 
prove some mathematical facts about the properties of the 
state transition matrix. 

The first small proof is quite simple: 


0 ( k , k ) - exp ( 


k 

F ( t ) dt ) 

. t- k 


(VI .9) 


In words, a state transition matrix not proceeding in the 
forward or backward direction in time (and thus remaining 
at the instant of concern, t^) is equal to the identity 
matrix. The proof of a second small fact begins with state 
transition matrices going the opposite direction in time 
over the same time interval 


® ( k 1 , k ) 


0 ( k , k-t 1 ) 


exp ( 


t k + 1 

F(t. ) dt) 
. fc k 


exp ( 


ft k 

F ( t ) dt) 

J tkr 1 


exp ( 


tk+l 
tk 


F(t) 


dt ) . 


From the above equation, it is apparent that 


0 ( k , k + 1 ) = 0* 1 ( k h - 1 , k ) 


(VI . 10) 


Equation (VI. 10) illustrates the fact that for a physical 
system of given state transition matrix, the state transition 
matrix going backwards in time, 0 ( k , k • 1 ) is equal to the 
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inverse of the state transition matrix g o i n g forward in time 
over the same time interval. Relations (VI. 9) and (VI. 10) 
illustrate the properties of the s t a t e transition matrix, 
and may help the reader to become familiar with the dynamic 
system models used in Kalman filtering/smoo thing. 

Now that one has found the discrete state transition 
matrix 4>(k+l,k) from the continuous matrix F(t), one should 
consider the case of a dynamical system with random inputs 

x ( t. ) = F ( t ) x ( t ) •+ G ( t ) w ( t ) . 

Please observe that w(t) is a continuous variable differing 
from the discrete noise process, w ^ . The subsequent, deriva- 
tion uses the fact that 

d 

<D(t,t 0 ) --- F(t)*(t,t 0 ). (VI. 1.1) 

dt 

which can be found from equation (VI. 8); and also uses the 
property of the state transition matrix 

<Mt 2 ,t 0 ) = •(t 2 .t 1 )*(t 1 .t 0 ). (VI .12) 

Relationship (VI. 12) can be shown to be true in the method of 
Ge lb et al . ( 1974 ) 

x(t 2 )=®(t 2 .t 1 )x(t 1 )=®(t 2 .t 1 )®(t 1 ,t 0 )x(to). 

But x ( 1 2 ) = ® ( t 2 . t 0 )x ( t 0 ) , 
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so that 4» ( t 2 . t 0 ) = <M t 2 . t J )® ( t I , t 0 ) Q . E . 1) . 

Or verbally, the state transition matrix over multiple time 
intervals can be represented as the product of the state 
transition matrices of each interval (this does not include 
random inputs) . 

The original derivation with random forcing follows 
Gelb et al. (1974). Starting with the fact that 
x(t) = <Mt,t(j) x(t 0 ) is the discrete solution to the 
homogeneous differential equation (VI. 2) with F time 
dependent, one proposes that the solution to (VI. 1) is of 
the form 

x ( t ) - «D( t,t 0 ) |(t) (VI .13) 

One then substitutes (VI. 13) into (VI. 1) 
d 

{•(t,t 0 )g(t) }«P(t)*(t,t 0 )g(t)>G(t)N(t), (VI .14) 

d t 

and upon taking the derivative and using the relation (VI. 11) 
one finds 

F(t)#(t.t 0 )|(t)**(t.t 0 )i<t) 'F(t)4>(t,t 0 )!(t)+G(t)w(t) (VI. 15) 

The time indices t and t () (as opposed to t k and t k+3 ) are 
used here for convenience. Eliminating terms in (VI. 15) and 
invoking (VI. 10) one obtains the equation (as per Gelb 
et al . , 1974 ) 
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|(t. ) = 4»(tg,t)G(t)w(t), 


(VI . 16) 


which can b e integrated dire c 1 1 y to y i e I d 

't 

£(t) = £(t 0 ) + ®(to.t)6(x)w(x)dT, (VI. 17) 

J t 0 

Again following the method of Gelb et al . (1974), one 

inserts (VI. 17) into equation (VI. 13), and noting that 
§ ( t q ) - x ( t o ) , thereby arrives at 

't 

x ( t. ) = 0( t , t 0 )x( t 0 )+ 4>(t,to)<MtQ,t)G(t)w(T)dx, (VI. 18) 

J to 

which simplifies using relation (VI. 12) to 

't 

x ( t ) = 0( t , to ) X ( t 0 ) *- 0 ( t . x )G ( x ) w( x ) dx . (VI. 19) 

J t 0 

Most apparently, if w(x) is zero (no random forcing!), then 
equation (VI. 19) reduces to the homogeneous solution x(t) = 
®(t, t 0 ) x ( t o ) , as should be the case. 

The second term on the right-hand side of equation 
(VI. 19) represents the stochastic contribution to the state 
parameter x(t). Equation (VI. 19) can also be written for 
the time interval t k to t k + i (Gelb et al., 1974) 

x(t k + i )=0(t k+ i , t k )x(t k ) + (VI. 20) 

ft k + l 

® ( 1 k + 1 > T ) n ( x ) w ( x ) d x , 

hk 
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and if one uses the foJ lowing definition 


I’ (k ) w k 


*-k +■ 1 

4> ( t k t .t)G(t)w(t)dx, 

. t k 


one can write (VI. 20) as 


(VI . 21 ) 


x k+l " ®(k+l,k) x k + r k w k , (VI. 22) 

using slightly modified notation. Equation (VI. 22) is the 
discrete (difference equation) form of continuous equation 
(VI. 1), and is also equation (III. 5), which is what was to 
be derived in this section. The pair of equations (VI. 22) 
and (III. 4) are the models on which one builds a Kalman 
filter/smoother. In general, one starts filter design with 
knowledge of the continuous dynamics of a problem as given 
by equation (VI. 1) and then forms more useful discrete equa- 
tions like (VI. 22) to work in combination with discrete meas 
urement data . 

Before proceeding further with the polar motion problem, 
it is proper to describe the two noise covariances, R(k) and 
Q(k), which are needed in the structure of the Kalman filter 
R(k) is the leasurenent error covariance matrix, while Q ( k ) 
is the more inscrutable system (or process) noise covariance 
matrix. Recall from (III. 6) that 

E{w k w 1 T ) = Q ( k ) 6 kl (VI .23) 
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E{v kVl T > = R(k) 6 kl (VI . 24) 

arid R ( k ) > 0. (VI. 25) 

where 6 k j is a Kronecker delta. R(k) is the more easily 
explained covariance, so it will be discussed first. The 
covariances R(k) and Q(k) are used in the Kalman filter as 
opposed to the noises v k and w k . 

The measurement covariance R(k) is related to the uncer- 
tainties v k in the measurements y k . Equation (VI. 24) 
suggests how R(k) is calculated. In the case of VLBI delay 

measurements, the Kalman filter used herein handles one 

2 

measurement at a time. Thus R(k) = v k for the scalar case. 
Also notice that R(k) > 0 for any measurement, which means 
j that no measurement is perfect (i.e. has no uncertainty). 

The standard deviation in the delay measurements is 

I 

equal to v k and is given by (Ma, 1978; Whitney, 1974) 

1 

v k = . (VI. 26) 

Ci> s SNK 

where w s is the "... standard deviation of the observing 
frequencies used to sample over the observing bandwidth...", 
and SNR is the signal to noise ratio of the observation. The 
signal to noise ratio is found from (Ma, 1978) 

[ 

2P 1/2 
SNR = (2BT) ' , 

rc 

where B is the recorded bandwidth, T is the integration time 
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1/2 

anci p^=Y(K) According to Ma ( 197 8), Yls the fringe visi 

bility and is equal to one "... for a completely unresolved 
source . . . " and 


I'al T a2 

K = 

(T a i + Tsl)(Ta2 + T s2 ) 


where 

T a j , T a2 = source antenna temperatures 

T S1 , T s2 = system [receiver! noise temperatures. 

The 1 and 2 subscripts denote the antennas at either end of 
the VL8I baseline of interest. Furthermore, the source 
antenna temperatures are defined by (Ma, 1978) 


T 


a 


FEA 


2k 


where F 

E 
A 
k 


source flux density 
antenna efficiency 
[antenna] collecting area 
Boltzmann's constant 


The VLBI delay uncertainties are typically calculated for use 
in least-squares analysis and stored in VLBI observation (com 
puter) files, so they are taken directly from these locations 
for use in Kalman filtering. The measurement covariances 
R(k) used herein do not include the effect of unknown 
systematic errors which could introduce biases into the 
estimates . 


104 



The determination of Q(k) is not as straightforward as 
that for finding R(k). in fact, it is not just Q(k) 
that is required in the Kalman filtering step (III. 9), but 
rather I’ ( k ) Q ( k ) I ,T ( k ) . The quantity r ( k ) Q ( k ) I ,T ( k ) is 
found using equation (VI. 21) 


l ’k w k 


-t k+l 

4>(t k + 1 , x )G ( x ) w( x ) dx . 

J t-k 


( VI . 27 ) 


Taking the transpose of the above, and changing the dummy 
variable from t to a results in 


T T 
w k r k 


tk+l T T 

w («)G (a)® ( t k 4 . j , a ) da , 

. t k 


(VI . 28 ) 


and by multiplying (VI. 27) with (VI. 28) and taking the expec- 
tation (E) of each side of the product, one finds, with the 
aid of (VI .23) 


r k Q(k)r k T - E <i’ k w k w k r k T > 


( V I . 2 9 ) 


= E< 


t k-* 1 


Jtk 


' K 1 

( t k + j ,x)G(x)w(x)w r (a)G T (a)<t> T (t k + i ,a)dxda) . 

. t k 


The expectation operation is an integration process, 
and since the order in which the integrations are done in 
(VI. 29) should not. affect the ultimate result, one can write 
(VI. 29) as 
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r kQ ( k ) r k 


t k + l 
t k 


f k + 1 

0(t. kfl ,x)G(x) 

tk 


( V I . 3 0 ) 


x E { w ( t ) w r ( a ) } G T ( a ) T ( t k + i , a ) dt da 


The expectation of the continuous white noise product 
E{w(t) w T (a)} is well known (Gelb et al . , 1974) with the 

value of Q(t) 6(x - a), where 6(x - a) is a Dirac delta 

function. Placing Q(t) 6(x - a) into equation (VI. 30) 
and integrating over alpha, one concludes (Gelb et a 1 . , 
1974 ) 


r k q(k )r k T --- 


t k+ 1 

® ( t k + 1 

t k 


X ) G ( X ) Q ( X 


)G T (x)4> T (t k) . 1 


x ) dx 


(VI . 31 ) 


Equation (VI. 31) relates the system (process) noise covari- 
ance Q(k) to a (continuous) spectral density matrix, Q(x). 

For more details, see Gelb et al. (1974). The term 
T 

rkQ(k)r k is inserted directly into the Kalman filter 
equations. In the research performed herein l' k is set to 
identity, I (for convenience), and any proportionality 
constants are absorbed into Q(k), the system covariance 
matrix . 

While it is relatively easy to estimate the uncertainties 
in delay measurements v k , and thereby the measurement covari- 
ance R(k), it is not so easy to fix the values of the system 
noise w k and its associated covariance, Q(k). People are 
sometimes critical of Kalman filtering because they don't 
understand how Q(k) is determined. It is sufficient to say 
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that there are severaJ methods that can be used to find Q(k), 
but it should be understood that wj< and Q(k) are strictly 
related to the physics of the system being modelled, and as 
such, can be found at least, in approximate fashion. For a 
strictly deterministic system, Q(k) = 0. The topic of wj< 
and Q ( k ) determination will arise at several subsequent 
points in this work. An example should help highlight one 
method of how wj< can be arrived at. More exact methods are 
available. 

Let one consider, as is done in VLBI , how time is kept 

by a maser. Conventional hydrogen masers presently keep 

time to one part in 10 44 or 10 45 . This is approximately 

one nanosecond in one day! This number does not account 

for deterministic variations in maser time due to expansion 

of the maser oscillation cavity, electrical problems etc. 

It makes sense to believe that the stochastic variations in 

1 4 

maser time w^ are on the order of one part in 10 . Thus 

the system covariance matrix can be approximated crudely as 
Q(k) * ( 1 . 0/1 . OxlO 44 ) 2 . More sophisticated maser time 
("clock") models will be discussed later. It should be 
pointed out that for estimation of some parameters, such as 
the time kept by maser, the parameter estimates are only 
weakly dependent on the exact choice of system noise covari- 
ance Q(k), and thus approximate values of Q(k) can be used in 
the Kalman filter. Finally, one does not always know the 
exact system noise information to use in the filter. 

Now that the covariances have been described, the 


107 


details of the system equation ( i . e , equation I l I . 5 ) for the 
various parameters of interest in VLBI -polar motion research 
will be set forth. Much should be known to the reader 
already about the measurement equation (III. 4), as the 
delays y^, the uncertainties in the delays v^, and the meas- 
urement partials M(k) have been described earlier; as such, 
the measurement equation will not be particularized further. 
On to the system equation. The vector of parameters to be 
estimated (x^) using VLBI delay data in this work is 


x k 


t c ( 1R) 
t cr ( 1R) 

L z 1 

t c (2R) 

x cr (2K) 

L 7. 2 


(VI .32) 


x 

y 

UTl -TAI 
A^ 

l Ae 


where x, y and (JT1-TAI are the polar motion and rotation 
components of primary Interest, and ^clock/ratef^) are 
"clock" parameters for the relative time between VI.BJ station 
i and some reference station R. One VLBI observation site is 
chosen as the reference station in each VLBI data set. and 
A6 are adjustments to the nutation series. (L z j = L z j , L z2 . 
L Z 3 , ...} are the atmospheric (wet zenith) delays as esti- 
mated at each station j. The number of clock and atmosphere 
terms to be estimated depends on the number of antennas 
participating in the VLBI observations. Also, the position 
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(declination and right ascension) of the quasars and radio 
sources observed in the VI.IJ1 data runs and the terrestrial 
coordinates of the stations participating in the observing 
will not be estimated as a part of x ^ . Instead, these posi- 
tions will be taken from multiple data set, least squares 
analyses of VLBI data done previously; these solutions are 
updated periodically by the VLBI personnel at NASA GSFC. 
While this may not be the best moans of fixing the quasar 
and VLBI station positions, such a procedure is necessary 
because of computer limitations, and has been shown in the 
past to be a reliable way to arrive at accurate parameter 
estimates. 

The discrete, system dynamical equations will now be 
derived for the polar motion components, mj and m 2 . Recall 
from equations (II. 7) that simple Chandler wobble (with no a 
priori excitation function forcing; 4* 1 . 4* 2 = 0 ) is 
described by the coupled relations 

Pij = O ( VI . 33 ) 

m 2 - o m i , (VI. 34) 

where o r has been replaced by the actual Chandler wobble 
frequency, o. The excitation functions are being set to 
zero so that one is assuming no pre - conceived form of the 
excitation functions which are of geophysical interest. The 
VLBI astronomical measurements will be used to infer what 
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may be happening to drive the Chandler wobble of the earth. 
Equations (VI. 33) and (VI. 34) can be written in matrix 


form as 



/ 0 

i 

— l 

\ * 2 / 

\o 



j x, (VI .35) 


an equation which does not account for any stochastic vari- 
ations in n»i and m 2 . Relation (VI. 35) clearly resembles 
equation (VI. 1) 


x ( t ) - F ( t ) x ( t ) 


with no stochastic forcing terms. By comparison of equations 
(VI. 1) and (VI. 35), one can conclude that 


F(t) = F 



( VI . 36 ) 


and that F is time invariant. If F is time invariant, one 
can f ind the discrete state trans ition matrix <D(t k + i, t k ) 
from F by equation (VI. 5) and (VI. 6) 

<F(t k+1 -t k )> 2 

0(k+l ,k)=exp{F(t k+ i t k ) } - I+F(t k + J t k )+ + 

2 ! 


Inserting (VI. 36) into equation (VI. 6) one finds (with At s 
* k +■ 1 ' 1 k ) 
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( oAt ) ' 


( oAt ) 


( oAt ) 


3 


<t> ( k *• 1 , k ) = 


1 - — 


-oAt + 


2 ! 


(OAt) 


(oAt) 


3 ! 


4 ! 
3 


(oAt ) 


2 ! 


3 ! 
2 


(oAt ) 


4! / 


or in more familiar form 


® ( tk+ i . tk ) - 


( cos(oAt), -sin(oAt) 
sin(oAt), cos(oAt) 


(VI .37) 


The discrete form of the polar motion system dynamics with 
no stochastic inputs can thus be written 



( cos ( oAt ) , 
s i n ( oAt ) . 


-s in ( oAt ) 
cos ( oAt ) 


m i 1 
m 2/k 


( VI . 38 ) 


In VLB I work. At. is ori the order of minutes, thereby making 
the argument oAt small; thus 4>(tj t+ j, t^ ) Is approximately 
the identity matrix. The sinusoidal form of the state tran- 
sition matrix makes sense since the original continuous dy- 
namic relations (VI. 33) and (VI. 34) model a coupled oscil- 
lator. 

Equation (VI. 38) is incomplete and now requires stochas- 
tic input terms. The polar motion process noise is assumed 
to be a random walk (or integrated white noise) which can 
be described by the equation 


x ( t ) = G ( t ) w ( t ) 


( VI . 39 ) 


or more explicitly 
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X 


(VI . 40 ) 


/■l\ /I 0 

\ m 2 ) \° 1 

where Wj and w 2 are white noise. G(t) is obviously identity 
here. The covariance of this integrated white noise . process , 
Q'(k), can be found from equation (VI. 31): (Q 1 is used for 
notation to represent the process covariance with the factors 
r k I’ k T absorbed . ) 



Q- ( k ) =I' k Q ( k ) r k T 


ft k+ l T 

®(t k + 1 .x)G(t)Q(t)CI 1 (T)®(t k + 1> t)dt. 

J * k 


(VI . 41 ) 


The components of the system noise spectral density matrix 
are defined as 


Q(x) = 


911 

912 

921 

922 


(VI . 42 ) 


With the preceding definitions inserted into equation (VI. 41) 
and ®(t k+ i, t k ) as defined by equation (VI. 37), one finds 


Q’ (k) 


9n' Q 1 2 ' 
921 ' 922 ' 


( VI . 43 ) 


' t k + i /cosoAt, -sinoAt\ / q^j 9 l 2 \/ c08 °At, sinoAt\ 
t k \sinoAt. , cosoAt j lq 2 i 922 /\“ s * n °At • cosoAt/ 


dt , 


where At now is t k + 1 -t. One can approximate equation 
(VI. 43), with ® ( t k H . i , t ) » I, as 
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Q' (k) * 


(VI .44 ) 


t-k + 1 , 

flu 

Ql2 

* 

^ <121 

d 2 2 


dx . 


A further approximation made in the filter model for the 
polar motion stochastic (process) covariance is that the 
mj and m 2 spectral densities are not correlated, 
thereby giving 


Q' 00 


t k+ l 

/ qn 

0 

*k 

l 0 

Q22 / 


dx . 


(VI .45) 


If qn and q 22 are constant over the range of integration, 
one can further approximate (VI. 45) as 

/ ^ll < t k+ 1 _t k * 0 \ 

Q'(k) - (VI. 46) 

\ 0 *I22< tk + l-tk) / 

In general, if one is going to determine the numerical 
values of qjj and q 22 by a trial-and-error technique during 
simulations with VI.BI data, then the simple form of equation 
VI. 46 is conducive to this task. Eventually it will be 
demonstrated herein that the polar motion is quite deter- 
ministic (Q'(k) « 0) and thus the form of the covariance 
found here really doesn't matter. But it is always a good 
practice to set up a Kalman filter model with all parameters 
having non- zero process noise covariances. 

Now that <D(k + l,k) and Q'(k) have been determined for 
the polar motion parameters, the state transition matrix 
and process noise covariance will be posed for the earth 
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rotation parameter * 3 . From equation (11.8), again 
with zero excitation (4*3). one has the continuous state 
equation for rotation with no random inputs 

x = m 3 = 0 . (VI . 47 ) 

Comparison of (VI. 47) with equation (VI. 1) shows that 
P ( t ) = 0 and thus from relationship (VI. 6 ), it is clear that 
the state transition matrix for m 3 is 

« ( k 1 . k ) = I , (VI. 48) 

or one in the scalar case. Again, a random walk process 
is assumed for the stochastic contribution to the rotational 
component, m 3 (the same as UT1-TAI), and thus 


x = m 3 = w. 


(VI .49) 


The process covariance for the rotational component is 
simply 


Q’ (k) 


1 k+l 

Q(t)dx 


(VI . 50 ) 


J fc k 


since •<t k+ l.T) and G ( t ) are both identity. The assumption 
that Q(t) is invariant in time is also made and the resulting 
covariance is 


Q' (k) = Q(t k + 1 - t k ) . (VI .51 ) 

A stochastic equation of the form x = w (with «>(t k + j,t k ) ) 
= I) represents a parameter x whose value is deterministic- 
ally constant, but has random variations driving it, in this 
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case a random walk. According to Gelb et al. (1974), a 
"... random walk process results when uncorrel ated signals 
are integrated. It derives its name from the example of a 
man who takes fixed-length steps in arbitrary directions. 

In the limit, when the number of steps is large and the indi- 
vidual steps are short in length, the distance travelled in a 
particular direction resembles the random walk process.” 

There are many naturally occurring entities of interest to 
mankind that are randomly driven by integrated white noise, 
and the assumption of parameters being uncorrelated is 
common . 

One also has to estimate offsets to the nutation in 
obliquity and longitude with respect to the (reference 
value) IAU (1980) theory of nutation (Kaplan, 1981). 

Recently Herring et al. (1985) have discovered annual vari- 
ations of several marcsec amplitude from the 1980 IAU model, 
so these offsets are estimated using crude Kalman filter 
models here. The state transition matrix is set to identity 
and the process noise covariance matrix for nutation is 
similar iri form to equation (VI. 46) 

/ ‘111 < t k+ 1 ) - 0 \ 

Q ' ( k ) = I (VI. 52) 

\ 0 Q22< t k + l _t k) / 

Perhaps better models will be needed in the future, but the 
ones used here should suffice for the present in dealing with 
data sets of one day duration or less. Currently, there are 
insufficient data to make any reliable adjustment to the 
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precession constant. 


It is of interest to examine the system model for VLBI 
"clock” parameters which were detailed earlier. As stated 
previously, a VLBI clock is actually not just time as kept 
by a maser, but rather also bears the influences of the 
phase injection and cable calibration procedures. It may 
also be composed of other time like variations due to unknown 
physical phenomena and oddities in the VLBI observing system 
and experiment. It is a rare case when one can fully 
account for all the sources contributing to the clocks. 

In general a constant ramp function dominates the clock 
signals along a given baseline. This ramp function is due 
to offsets in frequency of the masers at each end of a 
baseline, and should be modelled deterministically. The 
clock behavior along a baseline is modelled using two param- 
eters, a clock offset x c (i.e. a bias) and a clock rate x cr 
(i.e. a ramp). The continuous system model describing the 
clock dynamics is as follows 





( VI . 53 ) 


and simply indicates that the derivative of the clock offset 
is the clock rate, and that the derivative (with respect to 
time) of the clock rate is zero. That is, the clock rate is 
reasonably constant in time as is true for most masers. The 
model could be made more general by allowing the clock rate 
to change in time, but such a model may be too elaborate for 
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present purposes, especially for the case of single-day VLB! 
experiments. 

Equation (VI. 53) can be readily made discrete using 
equations (VI. 5) and (VI. 8). The result is: 

/ x C \ /! t kfl -t k \ 

\ T cr/ k+1 \° 1 / 

This equation may seem overly simple to model all that has 
been discussed, but the stochastic terms have not yet been 
included. In addition, because the Kalman filter produces 
innovations (adjustments) to the clock parameters as it 
processes (passes through) the measurements, the optimal 
estimates of the clock parameters £ ( k + 1 | k + 1 ) will reflect 
the actual clock variations in the data. So, even though we 
are limited by our understanding of the "clock" system to 
simple models, the Kalman filter will account to some extent 
for the neglected variations. 

Now that the state transition matrix has been fixed, the 
stochastic fluctuations of a VLBI clock must be modelled. 
Recall that the clock being solved for in VLBI is actually 
the relative clock behavior at one VLBI station with respect 
to the clock at some reference station, and as such is pro- 
portional to the difference in the maser times at the two 
stations. The random variation in this time difference at 
the outputs of the hydrogen masers, on the time scales of 
interest in VLBI, takes the form of flicker noise (Allan, 
1983), which is a noise process with a power spectrum that 


1 c 
x cr/ k 


( VI . 54 ) 
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is inversely proportional to frequency. Power laws for, 
and illustrations of the form of, flicker and random walk 
noise as well as other fundamental noises, are displayed 
in Figure 11 (taken from Allan, 1983). It has been shown 
(Brown, 1984) that a flicker noise process can only be 
handled approximately in a Kalman filter, whereas the cases 
of white noise and a random walk can be treated exactly. 

This is due to the fact that the power spectrum S(f) of 
flicker noise is not inversely proportional to an even power 
of frequency. 

The clock covariances are modified from the work of 
Jones and Tryon (1982) and Tryon and Jones (1982) and take 
the form 


^'clock(k) 1 Qclock^k^-l l k) (VI. 55) 

Q ' clock rate(k) = ^clock rate(*k«-l"*k). (VI. 56) 

It was felt that these models should be satisfactory for 
describing maser behavior over a day since Jones and Tryon's 
study was for a maser intercomparison over several months. 

A random walk works well for "clocks" even though maser 
noise is flicker- like because maser timing variations are 
dominated by changes in signal path (in the VLBI electron- 
ics) rather than by noise in the maser. Thus, one really 
wants to model these environmental effects, and a random walk 
model for stochastic clock flue tat ions is probably as good 
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ORIGINAL PAGE IS 
OF POOR QUALITY 



their optimum predicted estimates. (From Allan 
1984 ) . 


as any other. It is certainly more physical than a white 
noise model alone. 

Very early on in the research presented here a comparison 
was made between least squares analysis of a VLBI data set 
and a crude Kalman filter run estimating a single clock pa- 
rameter only (no rate was estimated). In actuality what was 
being solved for was not purely a clock parameter but also 
included the influences of other parameters not being esti- 
mated, such as zenith atmosphere delays etc. The clock model 
employed an identity state transition matrix and an inte- 
grated white noise stochastic model. The results show that 
the residuals from each technique are comparable and of 
similar distribution and thus demonstrated that an early 
version of the Kalman filter was operational. The study is 
presented here for its instructional value. 

The example is a comparison of residuals on the baseline 
running between Algonquin Park (Ontario, Canada) and Fort 
Davis, Texas. The VLBI measurement session spanned two days 
and began on 24 August, 1984. The least squares delay re- 
siduals are shown in Figure 12a, while the Kalman filtered 
clock estimates are presented in Figure 12b. 

In Figure 12a, the time during which the VLBI mea- 
surements were made (in days) is displayed on the ordinate, 
while the delay residuals (nanoseconds) are on the abscissa. 
Each weighted VI.B1 datum is identified by a capital Roman 
letter. Data taken, but not used in the weighted least 
squares solution, are designated by a small Roman letter on 
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the figure. The letters are used to identify the radio 
source (such as 3C273B, 3C 345, etc.) observed by the VLBI 
antennas; the key to the right of the least squares plot 
defines the source associated with each letter. The plot 
shown has had clock offset and rate contributions removed 
from the delay residuals by least squares analysis. The 
resulting delay residual "fine structure" plot is sinusoidal 
in shape with an amplitude of 6 nanoseconds. The sinu- 
soidal shape may be artificial, and could be due to several 
"breaks" in the clock behavior. It is the VLBI analyst's 
job to fit polynomials to the sinusoid; these polynomials 
would approximate the VLBI clock behavior on the Algonquin- 
Fort Davis baseline. A very good fit to the clock curve 
would leave (as a remainder) a linear trend of delay resid- 
uals; the residuals would run (in time) along a line centered 
at 0.0 nanoseconds in the updated delay residual plot. The 
scatter of the newly fitted delay residuals should be 
considerably less than 6 nanoseconds, and typically is 
about 50 - 150 picoseconds. 

For comparison, the clock parameter estimates from the 
Kalman filter are shown as a function of experiment time 
(a running Julian date - in days) in Figure 12b. In this 
plot the quasar and radio sources are not designated by 
letter, but the clock offset estimates are denoted by 
crosses. The curves of Figures 12a and 12b are visibly alike 
allowing one to conclude that the Kalman filter estimates 
match residuals found by the least squares analysis. There 
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is a time offset in the resuits of Figure 12b which has been 


removed (by fitting) in Figure 12a. It would probably be 
more proper to compare the residuals of each technique, but 
residual plots from the Kalman filter were not available at 
this early itage of this work. Producing a Kalman filter 
which could solve for clock variations as a function of 
time (with limited analyst intervention) was a primary goal 
of this dissertation. 

The Kalman filtered clock curve of Figure 12b was arrived 
at using an integrated white noise covariance Q ’ ( k ) = 

Q ( t k + 1 " *k) as opposed to some flicker noise method. 

In addition, clock parameter estimates are relatively 
insensitive to the noise level magnitude of the random walk 
model, Q, as long as it is in the domain of physical reality 
for masers (and calibration delays) used in VLB! work. 

The clock contributions dominate the measured VLBI 
group delay residuals, in general, while the second largest 
VLBI nuisance factor arises from wet atmospheric path delays. 
The associated parameter being solved for using the Kalman 
filter is the adjustment to the atmospheric path delay in 
the zenith direction at each station. The dry part of the 
atmospheric path delay is usually not solved for but rather 
is adjusted using the Marini model dry atmosphere cali- 
bration which was mentioned earlier. 

The choice of the wet zenith delay deterministic model 
was limited by the fact that only one parameter for each 
VLBI observing station is set up in the SOLVE analysis 
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program for atmospheric delay estimation. Thus, the simple 
model of: 

<Mt k + i.t k ) - I (VI. 57) 

is used. In the future, one would anticipate estimation of 
an atmospheric delay rate term as VLBI data analysis becomes 
more sophiscated. For now it is assumed that the wet zenith 
delay adjustment is roughly constant, although again because 
of its Markov foundation, the Kalman filter can still track 
non-constant atmospheric changes, even if they are not 
modelled. The zenith wet delay covariance can be modelled 
by a random walk process, as has been shown by Davis et al. 

( 1985 ) . 

This concludes the derivations of state transition 
matrices and system (process) noise covariances needed to 
Kalman filter and smooth VLBI group delay data. The 
reader is probably wondering how all of the preceding fits 
together. This will be illustrated . First the measure- 
ment equation will be specifically written out (recall that 
the general form of the equation is y k =-• M(k)x k + v k ): 
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( VI . 58 ) 


( 8 x 9x 9x 9x 9x 9x 9x 9 x 
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Remember that x is the VLBI group delay and v(delay) is the 
uncertainty in the delay. The variables x, y, UT1-TAI , A'K 
A6 , x c , x cr and L z have been explained earlier 
(see equation VI. 32 and related text). Equation (VI. 58) 
simply relates the group delay t to parameters x at each 
observing epoch t^. The subscripts k have been left out of 
(VI. 58) to avoid clutter. 

The discrete system dynamical equation (i.e. equation 
111.5) is 
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(VI . 59) 



+ w k , 

where At = t k+ j - t k . Notice that if At is small, the 
state transition matrix, Q(k+l,k) for this problem approxi- 
mates an identity matrix. If one had no information about 
the physics or behavior of the VLBI parameters above, one 
could possibly start with an assumption of 0(k*l,k) = I for 
initial filter runs. Also, the fact that all the VLBI pa- 
rameters above are given non- zero stochastic noise terms 
(w / 0) is only to allow the filter form to be somewhat 
general. One may find that some of the VLBI parameters are 
purely deterministic. 

Equations (VI. 58) and (VI. 59) comprise the measurement/ 
system dynamics model employed herein for Kalman filtering/ 
smoothing VLBI group delay data. The previously derived 
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covariances, R(k) and Q(k), and measurement partials, M(k), 
in addition to 0(k+l ,k) from equation (VI. 5 9), are inserted 
in the Kalman filter/smoother algorithms. The net results 
of the filtering/smoot hi rig will be the optimal parameter 
estimates x and their associated error covariances P. 


128 


Chapter VII. 


TUNING THE KALMAN FILTER 

To tune a Kalman filter is to define the values of the 
stochastic noise covariances used in the filter. Generally, 
one starts with some initial idea of what the variability of 
a parameter is based on knowledge of the physical system. 

For instance, one could analyze water vapor radiometer data 
to discover the stochastic variation in wet atmospheric path 
delay along zenith. Once one has some inkling of the range 
of applicable stochastic noise, one can then use some other 
method to home in on more precise process noise values. 

The method employed herein is strictly one of trial and 
error. This would not be feasible to do if the number of 
parameters being estimated were large. The tuning of a 
Kalman filter is usually the most, inscrutable facet of its 
application. Some authors even characterize tuning as being 
more of an art than a science, although tuning always has a 
tie to reality in that the process noise level used in a 
filter must be physically believable. 

The ability to tune the Kalman filter used in this study 
by a brute force technique like trial and error is allowed 
by its limited realm of applicability. The VLB! data to be 
analyzed using the Kalman filter are the 1984-1985 IRIS 
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(International Radio Interferometric Survey) observation set. 


which consists of single-day VI. 13 1 measurement sessions run 
once every five days. The IRIS radio astronomy antennas are 
located at Westford, Massachusetts; Fort Davis, Texas (HRAS- 
Harvard Radio Astronomy Station); Richmond, Florida; 

W o 1. 1. z e 1 1 , Germany; and Onsala, Sweden. Since the Onsala 
observatory is not a full-time participant in IRIS, data 
from the Onsala baselines will not be analyzed. Also for 
the purposes of this study, the observing station at 
Westford will be the clock reference; that is, it will be the 
station whose masers other station masers will be compared 
to along any given baseline. This will simplify the tuning 
procedure . 

The trial- and-error tuning method has a very simple 
basis - one adjusts the process covariance values until 
the Kalman filter residuals y^l - M(k + 1) x(k+l|k) are 
minimized. For clarification, the preceding residuals are 
not the delay residuals y^ - y^. Most. Kalman filter 
tuning procedures are based on adjusting process covariances 
to minimize the Kalman filter residuals and their variances 
and to make the statistics of the residuals behave in a 
Gaussian manner. The residuals are examined along each 
baseline as opposed to all at. once in a group. 

The IRIS Kalman filter tuning was not done for all 
parameters at once, but rather was done for the parameters 
which have the greatest effect on the residuals first, then 
those parameters having the next greatest effect, etc., until 
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all tuning was completed. The parameters with greatest 
impact on VI.BI residuals are those of the clock offsets and 
rates, and the clock offset and rate process covariances 
along a given baseline were found simultaneously. The clock 
parameter estimates were made using data from the baselines 
Westf ord-Ri chmond , Westf ord-Fort Davis(HRAS) and Westford- 
Wettzell; relative clock behavior along any other baseline 
in the Westf ord~Ri chmond -HRAS-Wettzel 1 system can be found 
from combinations of clock solutions along the baselines 
which use Westford as a reference. What is specifically 
being determined by trial and error are Q(clock) and Q(clock 
rate) in equations (VJ.55) and (V1.56). Once these values 
are fixed, they will be used in the processing of the 1984- 
85 IRIS data. 

In order to find a Q(clock) and Q ( clock rate) that will 
work in processing all the IRIS data, one would like to 
tune by examining the residuals of a one -day data set which 
exhibits behavior typical of most of the IRIS observation 
dates. A data set $84MAR04XP (version 13) from the 
IRIS series was arbitrarily selected for the tuning opera- 
tion, mostly because it was the data set used in constructing 
much of the filter structure and software. The designation 
S84MAR04XP means that the data collection started on 4 March, 
1984 and that the data set was processed by the NASA Goddard 
Space Flight Center VLlll analysis group. The least squares 
analysis residuals of S84MAR04XP (ver. 13) seemed well 
behaved, which allowed one some confidence in this choice of 
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a tune-up data set. The tuning results are summarized in 
Table One. 

The tuning procedure along each baseline consisted of 
at least two steps: first, examine the average residuals 

(or their standard deviations) as a function of Q(clock), 
Q(rate) in four-decade intervals over the wide Q ranges 
Q(clock) = {10 3 1 , 10 1 ° seconds}, Q(rate) - {10 3 8 , 1 0 22 
seconds and choose a subsection of the Q ranges for 
further study: second, examine average residuals and/or 
their standard deviations over the subsection until 
{Q(clock), Q(rate)} are found to the nearest decade in the 
Q domains. Further specification to better than one power 
of ten of Q(clock), Q(rate) for each baseline is probably 
not warranted at this time, unless one would want to Kalman 
filter phase delay. In general, there are only insignifi- 
cant changes in the average residual or standard deviation 
over a decade in either Q value, in the region of minimum 
average residual. 

Figure 13 depicts the results of the coarse clock 
tuning procedure along the Westford-HRAS baseline. Contours 
of constant average value are drawn to highlight the 
structure in the plot. A region of minimal average 
residuals extends from the center of the graph to the upper 
right-hand corner along a line in which both Q(clock) and 
Q(rate) are increasing; the average delay residual varies 
little along this line. 

In Figure 14 a coarse plot of the "average residual" 
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Figure 14. Kalman filter standard deviations of 
residuals. Contour values are underlined. Plot 
is for tuning of clock covariance levels. 
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standard deviations is displayed. The structure is similar 
to the plot of average residual values, and contours of 
constant standard deviation are drawn. A region for further 
tuning is surrounded by a square made up of dotted lines. 

This area has been chosen for further tuning because the 
standard deviations of the residuals are close to minimum 
value and the residual standard deviations don't change 
much within the region surrounded by the dotted line. 

Similar plots are made for the finer (single-decade reso- 
lution) tuning study and ultimately the best Q(clock) and 
Q(rate) are chosen from such a graph. This tr ial-and-error 
procedure was also followed for the Westf ord-Ri chmond and 
Westf ord-Wettzel 1 clocks, and it was reassuring to find that 
the plots used to select Q values for each baseline all resem- 
bled one another. In addition, the Q values found in this 
manner have been applied herein to non-IRIS data with good 
residual fit results. 

The clock tuning done by trial and error yielded 

process noise standard deviations of about 10 100 pico- 

1 4 

seconds for clock offset and one part in 10 for clock 

1 4 

rate, which are both physically reasonable. A part in 10 
for the clock rate process noise makes sense because the 
VLBI masers are supposed to operate at this specification 
over the period of several hours. The standard deviation of 
10 - 100 picoseconds in clock offset agrees with current 
GSFC VLBI data analysis uncertainties in which the all- 
inclusive "clock" is a dominant factor. 
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After the process noise covariances for clock offset and 
rate were found, the process (or stochastic) noise level for 
the wet (atmospheric) zenith delays, a set of nuisance pa- 
rameters, had to be fixed by trial and error. Characteriza- 
tion of these delays at each observing location is important 
for VLB I data sets in which no (or bad) water vapor radio- 
meter (WVR) data are available, which at present, al most 
always seems to be the case. The implementation of WVR's is 
currently in its infancy. 

The random walk covariance level at each station ( HRAS , 
Wettzell, Westford, Richmond) was estimated using data from 
two individual baselines: HRAS -R i chmo nd and Westford- 

Wettzell. Along the Westf ord-Wet tze 1 1 baseline, the previ- 
ously determined clock models were used in solving for clock 
parameters in addition to wet zenith delay parameter adjust- 
ments at Westford and Wettzell. The wet atmosphere delay 
covariances were adjusted until the delay residuals were 
minimized. A similar procedure was followed for HRAS- 
Richmond. As with the clock tuning, the residuals were 

estimated using atmosphere delay covariance levels over a 

“18 -30 

coarse range of 10 to 10 seconds in four decade inter- 
vals. This range was chosen because it is on the order of 
the square of nanoseconds to picoseconds - the level at 
which one might expect to see atmospheric delay variations. 
Figure 15 depicts the delay residual standard deviations as 
a function of covariance level for each of the stations at 
the end of the HRAS-Richmond baseline. The dotted line on 
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Figure 15. Kalman filter standard deviations of 
residuals for atmospheric covariance level tuning. 
VLBI stations are at HRAS (Texas) and Richmond, 
Florida. Contour values are underlined. 
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the figure again indicates the region over which fine 
tuning (at one-decade intervals) took place, and from which 
the final covariance levels were found. 

The standard deviations of wet zenith group delay were 
found to be on the order of tens of picoseconds, which again 
agrees with the physical value of what one would expect as 
the wet delay uncertainty. It should be noted here that 
estimates of the wet zenith delay adjustment should be 
directly comparable in magnitude and shape to the delays 
output by water vapor radiometers. Now that the covariance 
levels have been set for the clock and atmosphere nuisance 
parameters one can now turn attention to the more physically 
important parameters of polar motion, earth rotation and 
nutation . 

The covariance levels for the two components of polar 
motion were found (simultaneously with clocks and atmos- 
pheres) using the method of trial and error for the VLBI 
triad Westford-Wettzel 1 -HRAS . The Richmond station was not 
used in order to save time during tuning -- the trial-and- 
error method becomes somewhat long-winded after a while. 
Residuals were only examined for the baseline HRAS-Wettzel 1 , 
which should be sensitive to polar motion. VLBI baselines 
trending North -South are generally sensitive to at least 
one polar motion component; for more about this topic, 
please refer to Lundqvist (1984). The polar motion covari- 
ance levels being determined are and q22 in equa- 

tion (VI. 46). 
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The results of the polar motion tuning are presented in 

- 1 8 

Figure 16. The covariance level range is 1.0 x 10 
-30 2 

to 1.0 x 10 radians /sec for each component of 
polar motion, and again the area chosen for further fine 
tuning is marked by the dotted line. The standard devia- 
tion of the group delay residuals is sometimes used in the 
process noise covariance estimation procedure rather than 
the average residuals because the spread in the residuals 
seems to be a more useful indicator in tuning than some 
average residual. The average residual is usually very 
close to zero. In addition the standard deviation of the 
average residual is positive in value. 

The polar motion process covariance levels are roughly 
the same for both components and very small in value. 

The small covariance levels indicate that the polar motion 
parameters are deterministic to a high degree. The non- 
zero polar motion process covariances found through trial 
and error are retained in the filter only so the filter does 
not "close up" on itself. A filter which closes up is one 
that no longer modifies parameter estimates even though new 
data are being processed. It is a condition in Kalman 
filters that should be avoided, even if by injection of 
slightly extra artificial noise. 

The tune-up of nutation and obliquity was similar to 
that of polar motion except that polar motion parameters 
were estimated also, in addition to nutation adjustments, 
clock terms and atmospheric parameters. The nutation in 
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Figure 16. Kalman filter standard deviations of 
residuals for polar motion covariance level tuning. 
Baseline is HRAS (Texas) to Wettzell (Germany). 
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obliquity and longitude adjustments also proved to be very 
close to deterministic. Better nutation models (determin- 
istic and stochastic) should probably be implemented in the 
near future. 

The last remaining covariance level to be found is for 
UT1-TAI, length of day adjustments. All previously men- 
tioned parameters were solved for and the residual standard 
deviation was observed as a function of the process covari- 
ance level, Q, in equation (VI. 51). The change in length of 
day is found to be very nearly deterministic. Thus all the 
parameters of interest (nutation, polar motion, earth rota- 
tion) appear to be quite deterministic while the nuisance 
parameters (clock, clock rate, wet atmosphere zenith delay) 
bear strong stochastic components. In the future, what is 
here modelled as stochastic may become deterministic as 
scientists become more familiar with the physics of clocks 
and atmospheres . 

Table One summarizes the process noise covariance 
information for the parameters estimated in this work. It 
is interesting to note that all atmospheric process covari- 
ance levels are the same except for the HRAS station at 
Fort Davis, Texas which appears to be less variable than 
its eastern U.S. counterparts. It seems that the general 
lack of water vapor in west Texas is exhibited in a low 
relative covariance level, at least for the day of data 
being used to tune the Kalman filter. In addition the 
Westf ord-HRAS clock rate covariance level is anomalously 
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Table One- Kalman Filter Covariance Levels 


COVARIANCE LEVELS (Westford = reference 


CLOCK LEVELS 


Baseline 

Westford-Wettzell 
Westford - Richmond 
Westford -HRAS (085 ) 


clock 

10 

10 

10 


offset 

4 § 6C) 

-23 

-23 


10 

10 

10 


clock rate 


30 

30 

32 


, -1 * 
(sec ) 


ATMOSPHERE DELAY LEVELS 

Station 

Westford 
We t tze 1 1 
HRAS (085) 

R i chmond 

POLAR MOTION 

Each component 

NUTATION 

Each component 

UT1 - UTC 

One component 


Level (sec) 


10 


-24 


10 


24 


10 


25 


10 


-24 


10 ( rad 2 /sec ) 


i 30 / ,2 , 

10 ( rad /sec ) 


_ o ft 

10 ( sec ) 


station) 
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(but not severely) low, suggesting that the relative maser 
behavior along the Fort Dav i s -Wes t f o r d baseline is better 
than along other baselines. 

It should be stressed that the process covariance levels 
of Table One should produce good -to- excel lent residuals 
if applied to 1984-85 IRIS data sets for which they were 
designed. However, as improvements in the Mark III VLBI 
data-taking systems are made in the future, one may have to 
modify the Kalman filter deterministic and process 
covariance models applied here. In add i t i on, f i ner tuning 
may be needed if one is to filter phase delay data rather 
than group delay observations. 

Apparently there is some precedent for this trial-and- 
error approach to Kalman filter tuning. Fitzgerald (1967) 
uses a trial and error Kalman filter tuning technique for 
aeronautical guidance problems. There are more sophis- 
ticated techniques for determining the process covariances Q 
and also for the measurement covariances R than the trial- 
and-error method. For details about these procedures see 
sections 8.8 through 8.11 of Jazwinski (1970). Simply put, 
these methods use predicted residuals r(k+l|k) and the 
constraint that the predicted residuals be consistent with 
their theoretical statistical properties, to estimate the 
process covariances used in a Kalman filter. A specific 
example (for 1=1) of a predicted residual is 

r ( k+ 1 | k ) = y k+1 - M(k+l)x(k+l (k) . (VII. 1) 

An advantage of this technique is that one is using obser- 
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vational data to tune the Kalman filter. Time considerations 


and computer limitations prevented implementation of such a 
filter for this work. The smoother requires no further 
tuning because it only operates on the output of the Kalman 
filter. 
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Chapter VIII. 


FILTER OPERATION AND ILLUSTRATIVE EXAMPLES 

This chapter explains the Manner in which the Kalman 
filter was used to process VI.BI group delay data. Some 
illustrative examples are presented to clarify what is going 
on. In general the filter can be set up to automatically 
generate the estimates, error covariances etc. which are 
used as input to the smoother. Once the filter run has 
finished, the smoother is simply turned on. Thus again the 
filter operation will be emphasized here. The smoother 
depends on a proper filter run in order to produce optimal 
results. 

So far, much has been described about the state transi- 
tion matrix 9(k+l,k) and the process noise covariance Q, but 
little has been said about the a priori estimates, Xg, and 
error covariances, Pg , used to initially start the filter; 
lack of emphasis on this subject is a distinct shortcoming 
in most Kalman filter texts. It is very desirable to 
start a Kalman filter run with the best a priori values 
possible. The source of such a priori's could be previous 
Kalman filter end point estimates which stop at the start 
epoch of the subsequent Kalman filter run, or could be 
measurements of a priori parameters made just before the 


145 


start of data taking. 

In parameter estimation of IRIS VI.BI data there are 
usually no prior measurement data since the preceding 
observing session was five days earlier. The end point 
estimates from the earlier session could be used as a prioris 
for the current filter run, but this is not done here 
because it is thought that the maser time offsets will 
change greatly in the (four) intervening days due to 
the clock ramp signals inherent in maser time keeping. 

Perhaps it is better to start a Kalman filter run with the 
belief that one has no information about the a prioris. 

This is the tactic employed in this research. 

To utilize the assumption that nothing is known about 
the a prioris, one sets xq (actually, the adjustments to 
parameters with respect to reference values) to zero and 
Pq to infinity. Instead of setting the error covariance, 

Pq , to infinity, it is sufficient to set Pq to 
some large value much larger than the parameter errors are 
thought to be. For the purpose of IRIS VLBI parameter esti- 
mation, Pq is assumed (for simplicity) to be a diagonal 
matrix with null values for off-diagonal elements, and with 
Pq diagonal elements given in Table Two. It can be seen 
that some of the initial covariance values are not overly 
large (i.e. approaching infinity), but the filter seems to 
quickly recover from this. 

Since the amount of VLBI data for single IRIS runs is 
limited to about one day duration, it was necessary to 
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Table Two 


Kalman Filter Initial Conditions 


Initial Conditions 


Parameter 

x 0 


Po 


Clock of f set 

0 . 0 

sec 

1 . OxlO -16 

2 

sec 

Clock rate 

0 . 0 


1 . OxlO -28 


Atmospheric delay 

0 . 0 

sec 

1 . OxlO -16 

2 

sec 

Polar motion 
(each component) 

0 . 0 

radian 

1ft 

30.0x10 

radian 

UT 1 - UTC 

0 . 0 

sec 

- 7 2 

5.0x10 sec 

Nutation 

(each component) 

0 . 0 

radian 

— 1 ft 

30.0x10 

radian 


147 


iteratively Kaimari filter the VLBI group delay data. Pq 
was set to the values listed in Table Two for all runs. 

For the first iteration, x 0 was set to zero. The ending 
estimates x from the first iteration were input as a prioris 
into the second Kalman filter iteration, and each iteration 
was repeated in a similar fashion. The iterative procedure 
was concluded when the estimated residuals showed minimum 
standard deviation. Usually this occurred after three or 
four iterations. One would believe that the number of 
iterations could be reduced by starting with larger Pq 
values, but there was no excess programming time to look 
into this possibility. In general, the standard deviation 
of the residuals found using the iterative Kalman filter- 
ing method matched those arrived at via least squares anal- 
ysis. This result validates the belief that the Kalman 
filter method can be used to eliminate analyst intervention 
necessary in VLBI optimal estimation using batch least 
squares. Kalman filtering on a "state of the art" computer, 
as opposed to the antiquated (Hewlett-Packard) HP -1000, 
would be a dream. 

Now that most of the preliminary details have been 
unravelled, it will be quite useful to familiarize the 
reader with some examples of Kalman filtering. The initial 
examples reflect the analysis of IRIS data set $85SEP20XI, 
version 11. Figures 17 a and b depict the clock (offset) 
and clock rate (filter) estimates as a function of experi- 
ment running time (horizontal scale) for the baseline 
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Figure 17. (a) Westf ord-Wettzel 1 clock offset 

as a function of time as produced by Kalman 
filter. Note convergence to actual clock offset 
at the beginning of the run. (b) Wettzell- 
Westford clock rate as output by the Kalman 
filter. Note that the relative clock rate is 
approximately, but not absolutely, constant as a 
function of time. 
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Wettzell to Westford. The clock offset plot Is clearly 

dominated by a ramp function which is actually the relative 

behavior of the Wettzell maser with respect to the Westford 

hydrogen maser; an offset in frequency of each maser gives 

rise to a ramp function in time, and the difference of two 

ramp functions is the ramp function shown. Notice at the 

beginning of the clock offset plot that it takes a short 

time for the Kalman filter to converge on the proper clock 

offset estimates. The clock rate portion of Figure 17 is 

not as continuously linear as the clock offset plot, 

but it can be seen that the clock rate is roughly constant 

(considering the small vertical scale of the plot) with 

- 1 4 - 1 4 

value between -206. x 10 to -220. x 10 seconds/second. 
The approximately constant clock rate estimate again 
indicates that the masers along the baseline Westford- 
Wettzell are working well. There are three distinct changes 
in slope visible in the clock rate plot which can probably 
be attributed to phase injection delays and/or cable cali- 
bration effects. The changes in slope on the clock rate 
plot suggest that the deterministic "clock" model (of the 
Kalman filter) which includes clock rates, but no clock 
acceleration parameter, is approximate. The model will work 
better in some situations than in others. 

Clock offset along Westford-Wettzell with much of the 
dominating clock offset and rate removed is depicted in 
Figure 18. That the clock rate has at least three distinct 
values is now clearly reflected in the finer scale clock 
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Figure 18 . Wet tze 1 1 -Wes t ford clock offset with 
primary linear clock trend removed. Variations 
in clock offset as a function of time can be 
observed . 
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offset estimates. This figure is included to remind one 
that fine structure can be hidden by dominant trends in 
data . 

The other "nuisance" parameters in VLBI optimal esti- 
mates are the wet zenith (atmospheric) delays (at Westford 
and Wettzell) which can be viewed in Figure 19. It is 
obvious from comparison of the "clock" plots and Figure 19 
that the clocks will affect VLBI observations more than the 
wet zenith delay. However, it should be remarked that the 
plots of Figure 19 are for a day when variations in wet 
zenith delay are quite small. The graphs of Figure 19 
give the moisture content adjustment in the zenith direction 
at each station from VLBI data alone; this illustrates the 
sensitivity of Very Long Baseline Interferometry as a 
measurement technique. 

Now that some examples about the nuisance parameters 
have been presented, one can examine the polar motion esti- 
mates displayed in Figure 20. The Kalman filter standard 
deviations (from the square root of the relevant error 
covariance elements) are shown as bars with the polar 
motion component estimates on the figures. The polar motion 
and UT1-TAI estimates are made with respect to reference 
polar motion and UT1-TAI values from BIH (Bureau de l'Heure) 
publications (BIH Circular D or Rapid Service documents). 
Since BIH data are only found at five-day intervals, the 
reference values are determined by fitting a curve to the 
BIH polar motion data at the epoch of observation and at 
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Figure 19. (a) Westford atmospheric (zenith) 

delay adjustment as a function of time (as 
output by Kalman filter). (b) Wettzell 
atmospheric delay. 
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Figure 20. (a) Adjustment of Y-component of 

polar motion as output by Kalman filter; 
associated errors included. Error are at the 
one mas level. (b) Adjustment of the X-compo- 
nent of polar motion and errors. 
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two BIH values on either side of this IRIS datum; one then 
interpolates to get polar motion reference estimates (used 
in the VI.BI program CALC) at each VLBI observing epoch. 

A look at the general features of Figures 20 a and b 
reveals that the adjustments in polar motion in time are 
very nearly constant, especially considering the fine scale 
of the plots. One can also readily observe the estimates 
converging as more and more data are processed, as should be 
the case. The polar motion error bars originally start with 
amplitudes of about 2.0 mi 1 1 iarcseconds which fall to values 
of 1.5 mi 1 1 i ar cseconds or less at the end of the Kalman 
filter run. These error bar values lead one to conclude 
that the Kalman filtered VLBI polar motion estimates are 
good to the one- mi 1 1 iarcsecond level at best. 

In Figure 20 b, one can distinguish a not so gradual 
change in the x-component (mj) of polar motion of about 
0.4 m i 1 1 i arcseconds from the time 0.3 to 0.4 days. There 
is no corresponding change in the y-component of polar 
motion and, of course, a change of 0.4 m i 1 1 i arcseconds in 
the x-component is not statistically significant. Neverthe- 
less, it is interesting to note that this change in the 
x-component is well correlated in time with a large after- 
shock of the Great Michoacan, Mexico earthquake of 19 Sep- 
tember, 1985. The aftershock event was of magnitude M s = 

2 7 o 

7.6, moment 2.4 x 10 dyne-cm with epicenter at 17.802 N 

latitude and 101.647° W longitude (depth = 31 km) and 

occurred on 21 September, 1985 at 1:37:13.4 UTC (Preliminary 
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Determination of Epicenters, USGS , September, 1985). The 
VLB I experiment started at 17:37 UTC on 20 September, 1985 
and ran for about one day, with the 0.3-0. 4 day time range 
(when the change in the x-component of polar motion 
occurred) lasting from 00:49 UTC to 3:13 UTC on the 21st of 
September. The time of the aftershock falls well within the 
time range 00:49 to 3:13 UTC. 

It should be remarked that variations in clock and atmo- 
sphere parameters did occur from 00:49 to 3:13 UTC, and that 
these variations could influence the polar motion x - compo- 
nent results. But if this is indeed the case, one would 
expect such clock and atmosphere contributions to appear in 
the y-component of polar motion also. Such behavior is not 
evident and suggests that the change in the x-component of 
polar motion may be real, even though not measurable to our 
statistical satisfaction. Future improvements in VLBI 
measurement precision should allow resolution of whether 
such an observed variation is real. 

A return to the original topic of presenting examples 
is warranted. The change in length- of -day estimates for 
data set $85SEP20XI . version 11 appear in Figure 21. The 
figure shows that UT1-TAI is quite constant over the one-day 
VLBI session, although a very slight slope is evident in the 
estimates once the filter has recovered from the UT1-TAI a 
priori values at a time of 0.07 days. The error bars of the 
estimates are quite large until a priori recovery, which is 
what one would expect. Following the recovery, the 
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Figure 21. Kalman filtered adjustments to UT1 - 
TAI (changes in length of day) as a function of 
time. Note decrease in error size as filter 
converges on optimal estimates. 
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estimates converge and the error in the estimates decreases 
to 0.1 milliseconds or less. The reduction in the UT1-TAI 
errors over one day is relatively larger than that observed 
for the polar motion and nutation parameters. 

The adjustments to the nutation in obliquity and longi- 
tude are revealed in Figure 22. The nutation offsets are 
estimated with respect to reference values given by the 
1980 IAU (Kaplan, 1981) nutation series. Again the nuta- 
tion in obliquity and longitude offsets are roughly 
constant, and converge as the filter processes more and more 
data. The nutation in obliquity estimates are determined to 
a much higher precision than those of the nutation in longi- 
tude; this behavior is commonly repeated in other Kalman 
filtered VLB1 data sets. A step observed in both nutation 
estimates at time equals 0.84 days is due to influences of 
"clock" parameters; it is evident that the clocks are re- 
sponsible because both nutation parameters are affected. 

The last items of educational interest are the residual 
plots for a final run of the Kalman filter, which are given 
in Figure 23. The residuals are broken up into groups by 
baseline and, on a given plot, they would ideally be 
expected to have zero mean and some statistically predict- 
able value of standard deviation. The residuals are dis- 
played versus the running time during the observing session 
(on the vertical scale) so a VLBI analyst can spot trends 
in the data. The residuals along the Westford-Wettzell 
baseline bear some remnants of the "clock" conduct shown 
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Fig. 22b. Time in days 


Figure 22. Kalman filtered nutation adjustments 
and errors. (a) Nutation in longitude (b) 
Nutation in obliquity. 
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Figure 23. Kalman filter group delay residual 
plots. Roman letters Indicate (quasar or radio) 
source being observed. (a) Westf ord-Wettzell 
baseline. (b) Westford-Richmond (c) Richmond- 
Wettzel 1 . 
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previously, but at a much reduced amplitude. In general, 
the clock parameters have successfully absorbed such clock- 
like trends . 

Little average bias can be observed in the group delay 
residuals of Westford-Wettzel 1 and Westford-Richmond , while 
a somewhat larger offset in average delay is visible for the 
Richmond-Wettzel 1 baseline. Perhaps the bias could be 
removed by a slight retuning of the Kalman filter, but in 
most analyses such time-consuming procedures were avoided. 

The standard deviations of the group delays are in the 120- 
picosecond range which is typical for the data handled in 
this work. No large amplitude trends are evident in the 
residuals of each baseline, emphasizing that the Kalman 
filter produces results consistent with those attained via 
other current VLBI analysis techniques. 

A few more illustrative examples involving data sets 
other than S85SEP20XI follow. Data from $84MAR04XP are 
used to compare filter estimates with Rauch-Tung-Striebel 
(1965) smoother results. Filter estimates for wet zenith 
delay at Westford and Wettzell are in Figures 24. The corre- 
sponding smoother results are brought out in Figure 25. 

The rough, discontinuous wet delay fluctuations quite 
visible in Figure 24 become much smoother and more contin- 
uous in Figure 25, as one might expect. The smoother 
does "smooth" the data in the traditional sense. While the 
filter and smoother estimates are not very different in 
amplitude, the smoother uses all the information available 
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Figure 24. (a) Kalman filtered atmospheric 

delay adjustments for Westford, Mass. (b) 
Kalman filtered atmospheric delay for Wettzell, 
Germany . 
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Figure 25. Kalman smoothed atmospheric delay 
adjustments for (a) Westford. Mass, (b) 
Wettzell, Germany. 
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in the VLBI group delay data, and should produce slightly 
more accurate estimates than the filter. 

A comparison is now presented of wet zenith path delay as 
determined by Water Vapor Radiometer (Figure 26) and Kalman 
smoothed wet zenith path delay (Figure 27) as estimated 
using VLBI group delay data. The comparison is intended to be 
more for instructional purposes than for strict quantitative 
use. The data set used in the example is $85DEC10X, version 
6, with observing antennas at Westford, Wettzell and Onsala. 
Since the Kalman filter/smoother was not tuned for use with 
the Onsala station, covariance levels were set to values 
which seemed to be reasonable, and residuals were examined 
to make sure they fit well. 

In Figure 26 one can see that the WVR data density is 
quite high, with a gap during the middle of the run, while 
the data density is less for the Kalman filtered/smoothed 
zenith delay results. Nevertheless, it is apparent that the 
WVR and Kalman smoother zenith delays are roughly comparable 
in form and amplitude. There is a bias of approximately 
2.0 centimeters delay between the WVR and Kalman smoother 
results which is probably due to the method by which the 
Water Vapor Radiometer measurements (i.e. brightness temper- 
atures) are converted to wet vapor zenith delays. However, 
some of the bias could be due to somewhat approximate filter 
tuning . 

The Kalman smoother estimates do not have a gap near the 
middle of the experiment as do the WVR results, which 
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Dec 10, 1985 


Figure 26. Water Vapor Radiometer (WVR) wet 
atmosphere path length as a function of time 
(at Onsala, Sweden; 10 Dec, 1985); (Courtesy of 
Goran Lundqvist). 
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Tiae in days 


Figure 27. Kalnan smoothed wet atmospheric path 
delay adjustments (Onsala, Sweden - 10 Dec. 1985) 
as a function of tine. 
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suggests that Kalman smoother wet zenith delays can be avail- 
able at times when WVR observations nay not be. Also one 
nay observe that the smoother and WVR plots "track” one 
another quite well during the initial half of the VLBI 
data-taking session, while the natch is not so good during 
the second half of the run. This can easily be explained 
by the fact that in the filter/smoother, the wet zenith 
delay is modelled to be deterministically constant as a 
function of tine (with no wet atmospheric rate model). At 
the start of the VLBI session, the smoother, even with such 
a limited model, keeps up with the wet zenith delay slope 
until this slope changes. At this point the f i 1 ter/snoother 
has to deal with a change in wet atmospheric delay rate (a 
zenith delay acceleration!) while it expects the wet zenith 
delay adjustment to be constant. The result is that the 
wet zenith delay estimates become poor while the filter/ 
smoother tries to adjust to the new wet zenith delay slope. 
Still the smoother does produce a descending limb of wet 
zenith group delays. 

The final example in this chapter is included to demon- 
strate that wet zenith delay adjustments can vary greatly 
during the single-day period of a VLBI observing session. 

The results can be viewed in Figures 28, and are from 
VLBI experiment $86APR06X . version 5. The stations 
involved in the experiment were at Hatcreek, California; 
Owen's Valley Radio Observatory (OVRO), California; and 
Mojave, California. The data are very dense because the 
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Figure 28. Kalman filtered atmospheric delay 
adjustments as a function of tine for (a) Hatcreek. 
California (b) Owen's Valley (OVRO), California. 
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purpose of this experiment was to allow VLBI scientists to 
use phase delay observations during analysis as opposed to 
group delays; the former are more precise than the latter but 
phase delay analysis is more difficult than group delay 
work . 

The OVKO and Hatcreek wet zenith delay estimates from 
the Kalman filter show some initial recovery from a priori 
values, then reach a roughly steady state zenith delay, and 
finally experience a turbulent transition to a second steady 
state zenith delay which is of larger amplitude than the 
first delay level. This sequence of events is due to 
movement of a storm front across California. Some of the 
wet zenith delay fluctuation observed during the turbulent 
transition between the (steady state) zenith delay levels 
is due to use of poor quasar/radio source positions in the 
solution during this tine interval. The wet zenith delay 
amplitude change during the onset of the storm front is 
about one nanosecond, which is quite large and should not be 
ignored for VLBI estimation purposes. The Kalman filter 
readily includes the change of wet zenith delay as a part 
of the solution, while current VLBI least squares data 
analysis usually models the wet zenith delay as being 
constant over an one-day VLBI experiment. 
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Chapter IX. 


PRINCIPAL RESULTS 

This research began with the hope that a large earthquake 
(■agnltude 8.0+) would occur during an IRIS and/or Crustal 
Dynamics Project VLBI experiment, and thus would allow one 
to look at the effect of such an earthquake on the earth's 
polar motion during a single-day IRIS run. Using the 
Kalman filter one could look at the polar motion estimates 
as a function of time and thereby examine the Influence of 
an earthquake during the one-day length of data. Some 
questions to be answered would be: Mere there any precur- 
sory phenomena to the event?; Is the effect of an earthquake 
dominated by the main shock?; Is the earthquake signature 
actually a kink in the polar motion estimates?; and so on. 

To the author's chagrin, however, no major earthquake 

coincided with a 1984-85 VLBI observing session. The 

largest event occurring during 1984-85 IRIS observing was 

the Mexican aftershock described earlier; and, as was stated 

previously, any effect of this earthquake on polar notion 

was smaller than the polar notion error estimates and 

(therefore) not significant. It is indeed fortunate only for 

the purposes of this study that two large earthquakes with 

2 8 

moments of about 10 dyne-cn happened at all during 
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1985. These were the Nlchoacan. Mexico earthquake of 19 
September. 1985 and the event near the coast of central 
Chile on 3 March. 1985. At the 1986 Spring neeting of the 
American Geophysical Union, Adaa Dziewonski of Harvard 
University stated that these Chilean and Mexican earth- 
quakes are the first events of aoaent greater than 1.0 x 
28 

10 dyne-ca since 1977. 

Since these large events are available for further study, 
and such events are few and far between, this work has been 
relegated to deteraining if the aajor Mexican and Chilean 
earthquakes have had any effect on polar notion using the 
available IRIS data with five-day saapling interval. For 
both the Mexican and Chilean events, ten IRIS data sets 
surrounding each shock (see Table 3) were Hainan filtered to 
deteraine the polar notion coaponents on each of the ten 
dates . 

The filtering nethod applied was the iterative approach 
highlighted in the previous chapter. The use of ten IRIS 
data sets gives a forty-five day window around each event, 
and allows an individual scientist the chance to process 
the data by hiaself -- the analysis of aore data night take 
too auch tine, and such heroics are left to the national 
VLBI analysis teaas. The forty-five day tine window night 
allow one to see breaks in polar notion as Mansinha and 
Saylie (1970) observed, especially since they found that 
breaks in polar notion typically occurred within ±20 days 
of an earthquake. But, one should be forewarned that the 
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Table Three- Data Sets Analyzed Using Kalian Filter 


Mexico data interval 

Data set Version 


S85AUG26X1 8 

$85AUG31XI 8 

$85SEP05XI 11 

$85SEP10XI 11 

$85SEP15XI 8 

$85SEP20XI 11 

$85SEP25XI 10 

$85SEP30X I 11 

$850CT05X1 11 

$850CT10XI 11 


Chile data interval 


Data set 


Vers ion 


$85FEB07XB 

13 

$85FEB1 2X1 

14 

S85FEB17XI 

11 

S85FEB22XJ 

11 

S85FEB27X I 

11 

$85MAK04X I 

11 

S85MAR14XI 

12 

S85MAR19XI 

11 

S85MAR24X I 

1 1 

S85MAR29XI 

11 
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Mexican and Chilean earthquakes of this effort are somewhat 
saaller than the Alaskan and Chilean earthquakes which took 
place between 1957.0 and 1968.0, the data interval which 
Mansinha and Snylie exaained. More coaplete studies of 
IRIS data should be done in the future when wore data are 
available to test further the hypothesis that large earth- 
quakes signicantly drive the Chandler wobble. 

The data sets processed in this study are listed in 
Table Three. For the Mexico results, 3714 group delay 
observations were filtered, while 3501 group delays were 
processed for the Chilean event. No saoothing of the 
results was done. The outcoae of filtering each data set 
is the adjustaent in both polar aotion coaponents at the 
ending epoch of each run. These adjustaents are Bade with 
respect to BIH reference values. 

The Kalaan filter (polar aotion) adjustaents are dis- 
played in Figures 29 through 32. These adjustaents are 
coapared with VLBI polar aotion adjustaent values as 
produced by least squares analysis. The least squares 
values were taken froa a BIH-VLBI selection option of the 
SOLVE analysis prograa froa which one can choose either BIH 
or VLBI polar notion results as reference values for further 
data analysis. The BIH reference was chosen for the work 
coapleted here and is given as the zero aaplltude (adjust- 
aent) horizontal axis on Figures 29 through 32. 

It is apparent froa perusal of Figures 29 through 32 
that the least squares and Kalaan filtered adjustaents 
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Adjustments to Polar 
Notion (X) (ias ) 
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0.00 


Tine in days 


Figure 29. Polar notion ( X-coaponent ) adjust- 
aents (Kalaan filtered and least squares esti- 
aates) as a function of tine. Tiae interval 
surrounds epoch of the 1985 great Mexican 
earthquake. The Kalaan filter and least squares 
adjustaents track one another with a slight bias 
BIH estiaates are at a level of zero nas on this 
graph. LSQ values circled. 
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Figure 30. Polar motion ( Y-component ) adjustments 
(Kalman filtered and least squares estimates) as 
a function of time. Time interval surrounds epoch 
of the 1985 great Mexican earthquake. 
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Adjustments to Polar 
Notion (X) (mas) 



Figure 31. Polar notion (X-conponent ) adjustments 
(Kalman filtered and least squares estimates) as 
a function of time. Time Interval surrounds epoch 
of the 1985 great Chilean earthquake. 
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Adjustaents to Polar 
Notion (Y) (aas) 



Figure 32. Polar motion ( Y-component) adjustaents 
(Kalaan filtered and least squares estiaates) as 
a function of tiae. Tiae interval surrounds epoch 
of the 1985 great Chilean earthquake. 
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track one another quite closely on the plot, as night be 
expected since they cone fron the sane data sets; only the 
data analysis differs. Biases between the least squares and 
Kalnan filtered estinates are also evident. In addition, 
considerable offsets of nore than several nilllarcseconds 
can be observed between the VLBI polar notion adjustnents 
and the B1H reference values. Such offsets are expected 
since the BIH incorporates data fron techniques other than 
VLBI, as well as VLBI, in forning their polar notion esti- 
nates. If error bars were indicated with the VLBI estinates, 
they would be of approxinate anplitudes of one nilllarc- 
second. The error bars are not placed on the figures to 
facilitate conparison of the adjustnents. 

Part of the bias between the Kalnan filtered and least 
squares polar notion estinates arises because the Kalnan 
filtered estinates are fron the end epoch of the IRIS run, 
whereas the least squares estinates are for the epoch 0:00 
UTC contained in the IRIS run which is usually about 
eighteen hours earlier. Thus in Figures 29 through 32, 
polar notion estinates at (slightly) tenporally offset 
epochs are being conpared, thereby introducing uncertainty 
into the interconparison. The uncertainty could be as nuch 
as 4 nas (Dicknan, 1987). The conparison is only intended 
to denonstrate that the Kalnan filtered polar notion esti- 
nates are reliable as was indicated by the tracking of least 
squares and Kalnan filtered results. The renainder of the 
bias could be attributed to the fact that nutation adjust- 
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nents were not made in the least squares analysis or that 
the Kalman filter is not perfectly tuned. Some of the 
tuning deficiency Bay be due to the fact that the VLBI 
■easureaent errors are usually larger than the theoretical 
formal errors. Differences in the way wet atmospheric 
zenith delays are treated in the tmo analysis techniques 
may also contribute to the biases. The average biases in 
the x and y components of polar notion for the Mexican 
earthquake results are 1.7 ± 0.6 marcsec and 1.8 + 1.2 
narcsec respectively, and 3.3 t 0.6 marcsec and 2.4 ± l.l 
marcsec respectively around the time of the Chilean event. 

The ultimate results of this dissertation are presented 
in Figures 33 and 34. Two figures are the climax of this 
extensive effort? This limited outcome might remind one of 
high energy physics research in which several years of work 
are required to set up an experiment which is run in only 
several weeks or months. In Figure 33 the Kalman filtered 
polar motion adjustments are added to BIH a priori polar 
motion values (see equation III. 19) to produce the total 
polar motion estimates in x and y. The total polar notion 
component estimates are plotted versus one another in the 
figures. Annual wobble or the secular trend of polar 
motion have not been removed from these estimates. 

Figure 33 depicts a polar motion curve around the tine 
of the Michoacan. Mexican event. The curve is quite 
circular in general. Straight line segments have been 
drawn between the polar motion estimate points. The polar 
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Figure 33. Polar motion (X 
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Figure 34. Polar motion (X and Y components) 
estimates as produced by Kalman filter. Estimate 
interval surrounds epoch of the 1985 great Chilean 
earthquake. The arrow indicates the tine of the 
main event. Polar notion formal errors are about 
one mas. 
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motion errors in x and y are typieaily 1 mas. Experiment 
dates are indicated at the end points of the curve, and an 
arrow indicates the approximate time of the great Mexican 
earthquake of 1985. A large pair of error bars of amplitude 
10 marcsec by 10 marcsec are drawn on the plot to illustrate 
the approximate noise level which was prevalent in the (BIH) 
data Mansinha and Smylie (1970) worked with. 

The most striking feature of Figure 33 is the kink in 
the polar motion curve which is sharply coincident with 
the time of the 1985 Michoacan, Mexico earthquake. The 
kink appears to be real and is of amplitude greater than the 
one-slgma errors (*1 mas), although it is only a little 
larger in size than corresponding three-sigma errors («3 mas). 
The time correlation between the main event and kink in 
polar motion is impressive. One does not even require the 
curve fitting procedure of Mansinha and Smylie (1970) to 
see that it exists. 

The kink in the polar motion curve of Figure 33 begins 
prior to the 1985 Mexican event and continues until after 
the earthquake. This would lead one to believe that 
large mass movements occur before, during and after the 
main shock and are responsible for the kink in the polar 
motion curve and change in polar motion path. One could 
find the amplitude of the pole shift by a curve fitting 
procedure. It is doubtful that the observed kink in polar 
motion is due to the elastic event alone. 

The polar motion curve of Figure 34 illustrates the 
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time period surrounding the 1 9 8 5 Chilean earthquake. 

There were a large number of foreshocks and aftershocks 
associated with the main event, many of which occurred in 
swarms. The curve of Figure 34 is fairly smooth with no 
kink evident at the approximate time of the main event 
(indicated with an arrow). There is no polar motion datum 
on 9 March 1985, and this hinders the ability to resolve a 
kink in the curve during the earthquake. It should be 
pointed out that the Chilean earthquake was of magnitude 
7.8 as compared with 8.1 for the Michoacan, Mexico event, 
and this may be one reason to believe that it is easier to 
observe a kink in the polar motion at the time of the 
Mexican e v e n t . 

One might notice that the ends of the curve in Figure 
34 deviate from a true circular (or elliptical) path. This 
has been emphasized on the plot by adding an approximate 
path (fit by eye to the interior data) at the ends of the 
curve. The actual polar motion path begins to diverge from 
circular form at times of about ±15 days from the Chilean 
main shock. These deviations may be manifestations of the 
effect of the Chilean earthquake sequence (and associated 
mass movements) on polar motion (i.e. kinks in the path). 
That the breaks occur within ±15 days of the Chilean main 
shock casts a favorable l ight on the conclusions of Mans inha 
and Smylie (1970), for most of their polar motion breaks 
were within ±20 days of a large earthquake. The author has 
noticed that the times of the Chilean kinks coincide fairly 
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well with the epochs at. which earthquake swarms (affiliated 
with the main event) begin to occur and subside (Preliminary 
Determination of Epicenters, U.S.G.S., P e b r o ary -March, 1985). 

IRIS polar motion estimates (source: IRIS Earth Orien- 

tation Bulletin A, No. 29, July, 1986) over the intervals 
of the 1985 Chilean and Mexican earthquakes are displayed 
in Figures 35 and 36 for comparison with Figures 33 and 34 
(Kalman filter results). Aside from small biases, the 
polar motion path in Figure 36 closely matches that of 
Figure 34. Figures 33 and 35 are also quite similar in 
shape except at the endpoints of 26 August, 1985 and 10 
October, 1985. A quick look at Figure 30 of the VLBI y- 
component adjustments will confirm that the end point least 
squares estimates are quite different from the rest of the 
Kalman filtered and least squares adjustments. This 
suggests that these (least squares) end point values should 
be questioned, since the rest of the group of adjustments are 
maintaining a relatively constant bias with respect to BIH 
reference values. In any event, Figure 33 appears to 
bring out the existence of a kink in the polar motion curve 
which is not so certain in Figure 35. 

The (Kalman filtered) polar motion results of Figures 
33 and 34 are deconvolution filtered using the method 
applied by Gross and Chao (1985; see also Backus and Gilbert, 
1970) to LAG EOS polar motion estimates. Deconvolving the 
VLB! data in this manner produces excitation functions, ^ ( t ) , 
in which one should be able to directly observe shifts in 
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Figure 35. Polar motion (X and Y components) 
estimates as produced by IRIS. Estimate interval 
surrounds epoch of the 1985 great Mexican 
earthquake . 
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Figure 36. Polar motion (X and Y components) 
estimates as produced by IRIS. Estimation 
interval surrounds epoch of the 1985 great 
Chilean earthquake. 


the position of the pole about which polar motion takes 
place. IT IS STRESSED HERE that the excitation results 
presented subsequently include polar motion data for the 
first ten or eleven data points in each time series analyzed; 
the remaining points in each excitation series are zero 
padding (i.e. zeroes added to allow fast Fourier transforms 
to be used in the deconvolution). These points are not based 
on VLSI data and thus should be ignored in examining the 
results. They are included here for completeness in data 
presentation. 

Prior to deconvolution, a least squares fit of an offset, 
trend and sinusoid was made to remove polar motion long 
term behavior from the optimal estimates, thus leaving data 
which only bears the effects of earthquake or other short 
time excitation. A single sinusoid was used to represent 
the superposed annual wobble and Chandler wobble paths 
because the data being fit is of short duration. The fits 
obtained using an offset, trend and a sinusoid for each 
of the wobble components were as good as those obtained 
using an offset, trend and two sinusoids. 

The deconvolution filter applied to the polar motion 
data (after fitting and Fourier transformation to the 
f requency domain) is (Gross and Chao, 1985) 

G*(co) 

V(co) - 5 , (IX. 1) 

j ! G ( co) I j + (aN ) 

where the asterisk indicates complex conjugation, and the 



bars denote the taking of a norm. aN wili be discussed 
shortly. G ( w) Is the earth's transfer function (Gross 
and Chao , 1985) 


°0 

G(w) , (IX. 2) 

°0 “ ^ 

where 


o o - (2 jc/T) ( 1 + i/2Q) . ( IX . 3) 

T is the Chandler wobble period, Q is the quality factor of 
the Chandler wobble and Oq is the angular frequency (Gross 
and Chao, 1985). The Chandler wobble period T was taken as 
434 days; Q was taken as 106. 

The equation relating the Fourier transformed polar 
motion M(w) and excitation function 4* ( w ) is (Gross and Chao, 
1985 ) 


M(w)=4>(w)G(w) + N(w) (IX. 4) 

According to Gross and Chao (1985) the N(w) term represents 
"... the frequency content of the noise in the measurements 
of m(t). (polar motion estimates]" If the N(cj) were missing 
from IX. 4, the equation would relate what is going on in 
frequency space of a convolution process. To find ^(«) 
one would simply divide 
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GJ) 


1 


M ( w) . 


G ( gj) 

and then Fourier transform to find ^(t)- in this case, 

the filter applied to the transformed polar motion data is 
1/G(w) = G*/| j G | | 2 . 

One can observe from equation (IX. 1) that the filter 
actually applied to the transformed polar motion data M(w) 
is somewhat more general. Notice that if aN is set to 
zero, the general deconvolution filter becomes the simple 
deconvolution filter. Also if aN becomes very large, V(gj) 
approaches zero. 

a N is an adjustable parameter which controls how much 
measurement noise is allowed to remain in the excitation 
function estimates. For aN small, much noise exists in the 
excitation functions but temporal resolution is high. For 
aN large, resolution is lost but the noise level is 
reduced. The trade-off between noise reduction and 
achievement, of fair resolution can both be reached by 
choosing some intermediate value for aN in the deconvolution 
filter. Once the filter has been applied to M(u>), the 
results are Fourier transformed to find the excitation 
functions in time, ^(t). 

The polar motion data of Figures 33 and 34 are handled 
in the aforementioned manner, and the results are displayed 
in Figures 37 and 38. The x and y components of excitation, 
and 4*y, are calculated for the Mexican and Chilean 
data using aN values of 0.22, 0.5 and 1.0 following the 


190 



Excitation Function 
( x-mas ) 



Figure 37. Excitation Function ( X- component ) 
versus time. Time interval surrounds epoch of 
the 1985 great Mexican earthquake. Excitation 
function errors are about 0.5 mas. (a) orN = 0.22 
( N is a deconvolution parameter) (b) aN = 0.5 
(c) aN * 1.0. Note changes in scale on the 
ordinate . 
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Figure 37. Excitation Function (Y component) 
versus time. Time interval surrounds epoch of 
the 1985 great Mexican earthquake. Excitation 
function errors are about 0.5 mas. (d) ot N = 0.22 
(e) <*N = 0.5 (f ) ON = 1.0. 
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Figure 38. Excitation Function ( X -component ) 
versus time. Time interval surrounds epoch of 
the 1985 great Chilean earthquake. Excitation 
function errors are about 0.5 mas. (a) oN * 0.22 
(b) aN = 0.5 (this plot has been detrended). 

( c ) a N =■ 1.0. 
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Figure 38. Excitation Function ( Y-component ) 
versus time. Time interval surrounds epoch of 
the 1985 great Chilean earthquake. Excitation 
function errors are about 0.5 mas. (d) orN = 0.22 
(e) a N = 0.5 (f) aN= 1.0. 
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Fig. 38f . 
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method of Gross and Chao (1985). The Fourier transforms 
were typically done with 16 data values. 

It can be observed in the figures that vertical 
exaggeration (and thus possibly measurement noise content) 
is large for aN small (aN = 0.22) as was expected. The 
data from the figures for aN=0.5 will be analyzed further 
here since 0.5 is the value chosen as realistic by Gross 
and Chao (1985). Please note that all subsequent 
conclusions are based on this assumption. Also be aware 
that a linear trend has been removed from i^x^) for 
the Chile data. 

The excitation data for Mexico (aN=0.5) show small 
changes in , 4*y which are more gradual in 
structure than a step function. Variations in ^ start 
before the main seismic shock and continue on after the 
shock. The rates of change in ^x’ 4* y seem to be greatest 
within ±5 days of the main Mexican shock. These changes in 
ip are at the sub -mi 1 1 iarcsecond level and are thereby not 
significant when compared with the formal error level of 
about 0.5 mi 1 1 iarcsecond in each component of vp- 

The Chile excitation functions bear changes which are 
much larger than those of the Michoacan, Mexico event. 

These excitation steps occur within *25 days of the main 
shocks and a net change in excitation of about -2.4 
mi 1 1 i arcseconds is found in the x-component of ip. This 
ip x change takes place in two steps. 2.4 marcseconds is 
significant when compared with a formal error of 0.5 milli- 
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arcsecond . 


Zheng Ying (1986) presents strain results before, during 
and after large earthquakes. He proposes that strain steps 
can be observed a few months before a large earthquake and 
one month after such an earthquake. These observations 
confirm to a fair degree what is found herein using VLBI 
data near the epoch of the 1985 great Chilean earthquake. 

Since the VLBI polar motion results were deconvolved to 
produce excitation functions, it is appropriate to also 
deconvolve polar motion error estimates for comparison with 
these results. This will only be done for the 1985 Chile 
polar motion errors. The results are presented in Table 
Four. The deconvolution changes average polar motion 
errors of 1.2 and 1.2 mas (in x and y respectively) to 
average excitation errors of 0.5 and 0.5 mas. The exci- 
tation shifts observed (around the time of the 1985 
Chilean earthquake) are significant, at least at the one- 
sigma formal error level. 

Theoretically predicted shifts and directions (in 
degrees East longitude) calculated by Chao and Gross (1987) 
for earthquake dislocations only are compared with VLBI and 
LAGEOS observations of excitation shifts and directions in 
Table Five. The Sumba data are from LAGEOS observations, 
while the VLBI results are from this study. Observation is 
not strictly comparable with theory in this instance since 
the theory only models the effects of the elastic event on 
polar motion. The ratio of observed excitation to 
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Table Four- Chile Polar Motion and Excitation Function 


Measurement 
Date (1985) 

Polar 

X 

Motion 

y 

Excitation 
x y 

Feb 

07 

1 . 05 

1 . 03 

0 . 21 

0 . 77 

Feb 

12 

1.11 

1 . 16 

0 . 26 

0 . 69 

Feb 

17 

1.13 

1.11 

0 . 30 

0 . 63 

Feb 

22 

1 .07 

1 . 05 

0 . 34 

0 . 56 

Feb 

27 

1 . 20 

1 . 26 

0 . 39 

0.52 

Mar 

04 

1 . 53 

1 . 48 

0 . 48 

0 . 42 

Mar 

09 

- 

~ 

~ 

- 

Mar 

14 

1 . 24 

1 . 18 

0 . 57 

0 . 34 

Mar 

19 

1 . 32 

1 . 09 

0 . 64 

0 . 26 

Mar 

24 

0 . 99 

0 . 99 

0 . 69 

0 . 22 

Mar 

29 

1 . 07 

1.11 

0 . 76 

0 . 19 
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Table Five- 

Comparison 

of Observed 

and Theoretical 

Excitations 

Observed values >£ x 


1^1 (marcsec) Arg(^r) 

Mexico 

-0 . 2 

0 . 3 

0 . 36 

-56.3° E 
long. 

Chile 

-2 . 4 

-0 . 3 

2.42 

7.1° E 
long . 

♦ 

Sumba 

- 

- 

40 . 0 

-54 .0° E 
long . 

Theoret i ca 1 

values 

| >jr| ( marcsec ) 

Arg(tfr) 


Mexico 


. 084 

-83° E 


Chile 


0 . 18 

110° E 


Sumba 


0.21 

160° E 



Ratios ( ^-observed/ ^-theory) 

Mexico 4.29 

Chile 13.44 

Sumba 190.48 

* - Sumba event data from Gross and Chao (1985), and Chao 
and Gross (1987) 
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theoretical excitation is calculated also. Results are for 


the 1977 Sumba earthquake and for the 1985 Chile and Mexico 
earthquakes . 

The observed excitation shifts and direction due to 
the 1985 Mexico earthquake are roughly in agreement with 
theory. The excitation directions are only rotated 27 
degrees from one another and | i|>| ' s are close in value for 
theory and experiment. This is especially true when one 
considers that the formal error in either component of 
(observed) excitation is less than one mas. One might 
characterize the Mexico event as one in which polar motion 
theory and experiment are roughly in accordance. 

For the Chilean event the observed and theoretical 
excitation directions and magnitudes diverge more strongly 
than for Mexico. The theoretical and observed excitation 
directions are about 103 degrees apart, and the calculated 
amplitude ratio is about 13.4. This suggests that more 
may be taking place geophysically during the great Chilean 
earthquake of 1985 than is predicted theoretically. 

The Sumba (1977) results resemble those of Chile more 
so than those of Mexico. The observed and theoretical 
excitation directions are 146 degrees apart (almost 
oppositely directed from one another). The amplitude ratio 
is approximately 190. As one might recall from earlier on 
in this work, Gross and Chao (1985) suggest that the extra 
(and differently directed) excitation may involve the 
decoupling (due to slab pull) of a slab segment from the 
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subducted oceanic lithosphere during the Sumba event, and 
may be due to subsequent movement of the broken off slab. 
Independent evidence (tectonic force calculations) that this 
" s 1 ab-pu 1 1 -and -break " event did occur near the time of the 
Sumba earthquake is given by Spence (1986). The similarity 
of excitation results for Chile and Sumba leads one to 
believe that a s 1 ab- pul 1 -and-break event has occurred in 
Chile also. If so, such slab movements associated with 
great earthquakes may contribute significantly to the 
Chandler wobble excitation budget. In addition, unsteady 
movements of subducting lithosphere may excite polar 
motion. 

Several other earthquakes potentially capable of 
exciting the Chandler wobble of the earth have occurred in 
1985-86: Adak, Alaska; Taiwan; and Kermadec Islands. The 
Kermadec Islands earthquake is of particular interest since 
it was a magnitude M s =8.2 event (EOS, 2 December, 1986), 
and it compares in size with the Mexican and Chilean great 
earthquakes. Unfortunately, research into the effects of 
these events on polar motion was not performed due to lack 
of resources and facilities. In a related topic, one 
could check if a decoupled slab of lithosphere was formed 
during the 1985 Chilean earthquake by making calculations 
like those of Spence (1986). This would provide some 
independent confirmation that s lab-pul 1 -and-break took 
place. 

Since the author was unable to Kalman filter VLBI data 
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around the times (Sept. -Nov., 1986) of the Kermadec Islands 
earthquake (20 Oct., 1986, M s =8.2: EOS, 2 December, 1986) 
and of the Taiwan earthquake (14 Nov., 1986, M s =7.8; EOS, 

13 January, 1987), data from the IRIS Earth Orientation 
Bulletin A (January, 1987) will be examined. These data 
are polar motion results produced by the National Geodetic 
Survey (NGS), and they were estimated using least squares 
analysis. The raw data can be observed in Figure 39 
(Pole position 1983-1986). The data of interest are shown 
near the end of the pole position curve (marked by the 
Roman letter B). The pole position estimates displayed in 
Figure 39 includes annual wobble and Chandler wobble 
components . 

The data of interest from Figure 39 are deconvolved 
using the method of Gross and Chao (1985), and the exci- 
tation function estimates are shown in Figure 40. A 
linear trend (-.062 mas/day) has been removed from the 
component results. The plot shows little net change 
in the excitation as a function of time. The detrended 4*x 
component shows a net (maximum) change in excitation of 3.8 
mas; this is very significant when compared to a typical 
excitation function error of about 0.5 mi 1 1 iarcseconds . 

The x-component of excitation starts to change directly 
after the main Taiwan shock, and increases greatly there- 
after. Prior to the main Taiwan shock, the x-component of 
excitation is quite constant. The correlation of the onset 
of x excitation change with the time of the Taiwanese 
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Axis (msec arc) 



Fig. 39. Pole Position from Nov 17, 1983 (A) 

to Nov 30, 1986 (B) at 5-day intervals. 

Circled pts: X,Y = SLR [Satellite Laser Ranging] 
value; open circles: no value. (From IRIS Earth 
Orientation Bulletin A, January, 1987). 
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Figure 40. Excitation Function versus time. 
Time interval is near the epochs of the 1986 
great Tonga-Kermadec and Taiwan earthquakes. 
Excitation function errors are about 0.5 mas. 

N = 0.5. (a) X-component (detrended) (b) Y- 

component . 
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event is striking. The Taiwanese main shock precedes the 
change in polar motion excitation, thus inferring that 
earthquake (and post -seismic mass movement induced) exci- 
tation of the Chandler wobble is taking place. The source 
of such post-seismic excitation might be due to movements 
of lithospheric slab and/or "...viscous relaxation of the 
underlying asthenosphere . . . " around the earthquake zone (Au, 
1980; Nur and Marko, 1974; Cohen, 1979; Spence and Turcotte, 
1979). In either case, mass movements associated with the 
main earthquake shock seem to be affecting the Chandler 
wobble. 

Finally some comment should be made about how much 
aseismic events (and seismic events) associated with the 
Mexican, Chilean, Taiwan and Sumba earthquakes contributed 
to polar motion excitation. The sum of the excitations 
associated with these four events (occurring from 1977-1986) 
has a maximum value of about 46.6 marcsec. If one then adds 
in the theoretically determined excitations (Smith, 1977) 
associated with the Chilean and Alaskan earthquakes of 1960 
and 1964, one arrives at a total excitation (related to 
only six great earthquakes) of about 103 marcsec. If the 
Sumba results are not included, the excitation level would 
be 63 marcsec. This excitation level is considerable; 
examination of an excitation data set of duration similar 
in length to the Chandler wobble damping time would allow 
us to draw more valid conclusions about earthquake- 
associated excitation of polar motion. In addition, the 
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level of polar motion excitation due to great earthquakes 
and associated mass movements may vary with time, and other 
sources of excitation could dominate solid earth mass 
effects at various times during earth history. 

The polar motion results derived in this work can be 
interpreted in a plate tec t on i c/geophy s i ca 1 framework. 

(This is not to say that atmospheric and oceanic phenomena 
are totally unrelated to Chandler wobble excitation). There 
is evidence that the Chandler wobble was excited following 
the 1985 Chilean and 1986 Taiwan great earthquakes. This 
excitation could be attributed to one or more geophysical 
phenomena: post-seismic movement of a lithospheric slab 

(Gross, 1985) and/or (post-seismic) viscous relaxation of the 
as thenosphere under the earthquake zone (Au, 1980). 

For the 1985 Chilean earthquake, there is also (VLB1) 
evidence of pre-seismic (before the main shock) Chandler 
wobble excitation. Thus, VLBI could be used as a (temporal 
and spatial) earthquake prediction tool. If VLBI data 
could be shipped and processed in a timely fashion, perhaps 
one could make predictions of seismic main shocks in a 
range of several days to about one month prior to an 
event. The pre-seismic changes in wobble seen in the VLBI 
data would also demonstrate the validity of Kanamori's 
(1977) work. 

A possible source of (pre-seismic) Chandler wobble 
excitation is the bulging of the lithosphere at a Chilean- 
like subduction zone prior to a main shock. Such bulging 
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could change the inertia tensor of the earth, thereby 
exciting the Chandler wobble. Calculations should be made 
to determine the theoretical excitation of such bulging. 
Figures 41 and 42 (Turcotte et al., 1978; Uyeda, 1979) show 
such (bending and) bulging which may precede great 
earthquakes such as the 1985 Chilean event. The amount of 
bulging may increase greatly just prior (several months 
before) the main shock (thereby driving a great thrust 
earthquake). Such bulging might not occur at a plate 
margin such as the one near Taiwan (refer to Figure 42) and 
thus explain why Chandler wobble excitation steps are not 
observed prior to the 1986 Taiwan earthquake. Pre-seismic 
excitation may also be due to ase.ism.ic fault slippage. 

The main emphasis of the preceding lines of thought is that 
most of the earthquakes of our planet occur at subductiori 
zones, and it appears that a considerable portion of 
polar motion excitation may be driven by subduction- 
related phenomena. 

More VLB 1 data are needed during time periods surround- 
ing large/great earthquakes to determine the actual effects 
of such earthquakes on the Chandler wobble. Kanamori (1977) 
presents a listing of the great (shallow) earthquakes of 
1904-1976; there are roughly nine such earthquakes per 
decade over this time span. Seismic events, and related 
major mass movements, are a significant contributor to the 
to the Chandler wobble excitation budget. 
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Fig. 41. Bending of an elastic or elastic- 
perfectly plastic plate in response to applied 
bending moment M, vertical force Q and horizon- 
tal force S, with hydrostatic restoring force 
applied by the seawater above and assumed 
fluid mantle below. Redrawn with simplification 
from Turcot te et al., 1978. 
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CONCLUSIONS 

The primary goals of this work were first, to show 
that a Kalman filter can be used successfully to analyze 
Very Long Baseline Interferometry group delay data with 
emphasis on handling clock variability; and second, to shed 
some light on the problem as to what degree earthquakes 
excite the Chandler wobble of the earth. The former tech- 
nological objective has been completed to a reasonable 
extent while the latter scientific objective continues to 
elude total resolution as it has so often during the past 
century. Many more IRIS observing sessions will pass 
before the latter problem is adequately solved. 

The Kalman filter design put forth herein does a satis- 
factory job of estimating clock terms, wet atmospheric 
zenith delays, earth orientation and rotation parameters 
and nutation offsets. The clock models can handle all but 
very large changes in clock rate. Very large clock rate 
changes could be accounted for either by making manual 
adjustments during the data processing or by implementing 
a clock model which includes clock acceleration parameters. 
With any luck, most large clock rate changes will be arti- 
facts found only in VI. BI data sets of the past. 
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Wet atmospheric zenith delays estimated using the 
Kalman filter were found to generally match the zenith delays 
determined from water vapor radiometer observations, except 
for slight biases. This confirms work originally performed 
by Tom Herring and Jim Davis at the Harvard/Smithsonian 
Center for Astrophysics. 

Polar motion estimates from the Kalman filter match the 
values found via least squares analysis to modest precision, 
usually better than the agreement of BIH results and VLBI 
least squares estimates. The sources of biases seen 
between the least squares values and Kalman filter estimates 
have not been isolated due to time and funding constraints, 
but the causes are thought to be differences in analysis 
methods and assumptions. It could also be that the filter 
is producing more objective results than least squares 
methods, since Kalman filtering of VLBI delays requires 
minimal analyst intervention once the filter has been tuned. 
Nutation and change in length of day estimates were also 
made using the Kalman filter, but are not presented here 
because they were not parameters of interest to this 
study -- they were estimated in order to produce reliable 
polar motion results. 

The trial and error method of tuning the Kalman filter 
took much time and computer cycles but did produce good 
results. More recently discovered tuning methods (than 
trial and error) based on the adjustment of statistics of 
filter residuals to meet theoretical expectations should be 
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employed in future undertakings. 

It should be obvious from the last several chapters that 
the Rauch -Tung-Striebel smoother was only rarely applied to 
VLB 1 group delay data. ‘’Why?” would be an appropriate 
question. The Rauch-Tung-Striebel smoother requires the 
storage of all parameter error covariance matrices from the 
forward Kalman filter pass. A typical Kalman filter run 
could easily generate tens- of -thousands to hundreds-of- 
thousands of covariance elements which are difficult to 
store in memory. For this reason the smoother was only used 
in limited (i.e. small number of estimated parameters) 
applications. Filtering the data in a backward time direc- 
tion (back filtering) after forward filtering is probably a 
better idea than smoothing because storage required for 
back filtering is minimal. 

Some reference is now made to general problems encoun- 
tered in this research. The largest single cause of grief 
was having to use an Hewl ett-Packard- 1 000 computer. While 
the HP-1000 is an excellent machine for computer instruc- 
tion, it presents a nightmare for number crunching. The 
ability to load large programs is limited, and one is 
required to add significant amounts of software to use the 
extended memory (EMA) features, which is the default opera- 
tional mode on many modern computers. Unfortunately, most 
of the NASA VLBI data handling and analysis programs (CALC, 
SOLVE, Data Base Handler) are written for the HP-1000 and 
the code is not very transportable. 
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There were several other forays into the absurd which 


wasted much time during this research. Several weeks were 
spent implementing an algorithm called the epsilon 
technique (Gelb, 1974) in an attempt to make early versions 
of the Kalman filter operate properly. One should be warned 
that the epsilon technique is probably not needed in Kalman 
filter implementation on a modern computer, and at best only 
offers an additional adjustable parameter with which one can 
vary filter output. Efforts should be directed at properly 
setting up the state transition matrices and process noise 
noise models rather than invoking the epsilon technique. 

The other major sink of time and effort involved an attempt 
to use a flicker noise generator for the betterment of Kalman 
filter clock models. As has been alluded to earlier, the 
clock noise is probably dominated by changes in signal path 
rather than by a (flicker) maser noise component, so sever- 
al months were "spent" on this activity. 

Many approximations were made in the course of this 
research, and it would be prudent to re-emphasize a few of 
them. In the clock process noise covariance calculation, 
terms involving the time interval t^-t-i - t^ of the state 
transition matrix (see equation VI. 29) were not included in 
the calculation. This is because the clock model applied 
herein was modified from the work of maser experts (Jones 
and Tryon, 1982). In retrospect, such an approximation 
should not be made, but the author was giving in to simpli- 
fications which made filter set-up easier. 
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An approximation forced more by the SOLVE program 
structure than by anything else was that requiring the wet 
atmosphere zenith delay to be constant. It would be more 
accurate to allow for deterministic rate variations in the 
wet zenith delay, as are actually observed in WVR data. 

There is a tradeoff in having such extra (rate) parameters 
in a Kalman filter: the more accurate parameter model will 

necessitate additional filter tuning. Still, atmospheric 
rate parameters are part of a more precise VLBI future. 

It is obvious that much more IRIS data will be needed 
before any ultimate conclusions about the effects of 
earthquakes on the earth's polar motion can be drawn. 

This may take some time, as no very large events occurred 
between 1977 and 1985; this lack of activity over long time 
intervals might lead one to believe that earthquakes are not 
the only source of Chandler wobble excitation. Still, no 
conclusion is firm until more is understood about the 
damping of the Chandler wobble. 

One kink observed in the polar motion curve of 1985 
appears to be real and extremely coincident with the great 
Michoacan, Mexico earthquake. Polar motion kinks are also 
found within +15 days of the large 1985 Chilean event. The 
current IRIS one-sigma formal error level is about one 
milliarcsecond in each component of polar motion. This 
is a factor of at least ten better than the BIH results 
studied by Mansinha and Smylie (1970). Thus, while M a n s 1 n h a 
and Smylie were able to study the effect of very large 
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earthquakes (Alaska and Chile of the 1960's) on polar motion 
with low precision, this study has examined the effects of 
somewhat smaller events at significantly higher precision. 
One is forced to conclude that more and better data 
are needed. Another outcome of this study is that Kalman 
filters will be available to look at polar motion path fine 
structure (in time) should a very large earthquake take 
place during an IRIS observing session. 

Excitation functions determined from the optimal esti- 
mates of polar motion show statistically non-significant 
excitation of the Chandler wobble by the 1985 Michoacan, 
Mexico earthquake. The observed excitation at the time of 
the Mexican event coincides closely (in amplitude and 
direction) with theoretical values found. The static 
deformation field models do not incorporate the effects 
of aseismic slip or other mass movements. 

Significant excitation of the Chandler wobble (at least 
at the one-sigma level of measurement error) is observed 
(using VLB I ) near the time of the great Chilean earthquake 
of 1985. The magnitude of the excitation change is about 
2.4 mi 1 1 iarcseconds as compared with a formal error of 
about 0.5 mi 1 1 iarcsecond . Such an excitation is much larger 
than theoretical seismic models permit. The excitation is 
also (roughly) oppositely directed from theoretical vectors. 

The 1985 Chilean earthquake appears to affect the 
Chandler wobble in a manner similar to the great Sumba 
earthquake of 1977, although the effects are not as large. 
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It may be that lithospheric slabs are being broken of! of 
the subducted lithospheric plate during the main seismic 
shock of such a large earthquake. Thereby, the decoupled 
slabs are moving after such an event and change the inertia 
tensor of the earth. Independent tectonic force studies 
verify that such a slab-pull-and-break mechanism is 
realistic. We will need to know the degree of Chandler 
wobble excitation by such slab movements and the frequency 
of such events in order to discover to what degree this 
tectonic phenomenon drives the Chandler wobble. If the 
excitation step of the 1977 Sumba event is actually forty 
milliarcseconds, then decoupled lithospheric slabs drive 
the Chandler wobble at a large percentage of the total 
observed excitation level. 

The IRIS VLB1 results around the time of the 1986 
Taiwan earthquake show that the time of the Taiwan main 
event is strongly correlated with polar motion excitation. 

A change in excitation level directly follows the time of 
the 1986 great Taiwan event. The excitation is probably 
due to post-seismic movements of lithospheric slabs 
decoupled from lithospheric plates, or due to relaxation 
of the asthenosphere near the fault zone. The preceding 
mechanisms are likely candidates for post-seismic Chandler 
wobble excitation. One source of pre-seismic Chandler 
wobble excitation is posed to be (changes in the inertia 
tensor due to) lithospheric flexure at subduction zones 
such as the one near Chile. Aseismic fault creep is 
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another possibility. 


The future of Very Long Baseline Interferometry studies 
of the earth's polar motion is certainly bright, but every 
effort should be made to optimize IRIS results. Surely 
water vapor radiometers will aid in minimizing the effect 
of wet path delays on VLB! observations. More must be done 
to improve other calibrations also. Perhaps IRIS obser- 
vations could be made daily (instead of at 5-day intervals) 
in the next several years also. 

As a final remark emphasizing something which has been 
pointed out before, the IRIS network of Westf ord-Wettzel 1 - 
Fort Davls-Ri chmond- Onsala is far from the best observing 
configuration for earth orientation and rotation. A grid 
utilizing Westford (or Canada), Wettzell, a southern Argen- 
tinian station and an antenna in southern Africa would be 
better. Of course, efforts are being made in this direction 
and perhaps someday economic and political considerations 
will permit such an array to exist. 


225 


I 


APPENDIX I. 

KALMAN FILTER DERIVATION 

The discrete Kalman filter will be derived. The abbre- 
viated treatment presented herein follows Gelb et al. (1974) 
The equations which form the foundation of the Kalman filter 
are ( I I £ . 4 ) , ( I I I . 5 ) and ( I I I . 6) : 

Vk = M k x k + v k (A1 ' 1 > 

x k+l = ®k x k + ^k w k (AI.2) 

Efw^} = 0. ElVfc} = 0 ( A I . 3 ) 

E < w k w l T }=Qk 6 kl • v^l 1 >* R k 6 kl • 

Equation (AI.l) describes the measurement process, 
equation (AI.2) models the dynamics of the parameter vector 
Xfc , and equations (AI.3) give the assumed properties of 
the measurement noise v^ and process noise w^. 
y k is the vector containing the measurement information, 
and M|< relates the parameters to be estimated to this 
information. 0^ is the state transition matrix which 
models how x^ +1 is related to x^. and I’ ^ relates 
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the process noise to x k ^ . Q k and R k are the 
process noise and measurement noise covariances. Capital 
letters usually denote matrices while small letters are 
reserved for vectors . 

The first step of the Kalman filter algorithm to be 
derived here shows the relationship between optimal param- 
eter estimates in time. Beginning with equation (AI.2) 
and taking the expectation value of this equation, one 
finds 

E{x k + 1 > = ® k E{x k ) + r k E{w k } . 

Carets are used to denote optimal parameter estimates (also 
recall that E { w k ) =0 ) . The result of taking the expec- 
tation is 


*k+l = ®k*k- ( A 1 . 4 ) 

Equation (AI.4) demonstrates how optimal parameter estimates 
are projected in time. In the notation of Jazwinski (1970), 
equation (AI.4) is written as 

&(k + l|k) - 0» ( k+ 1 , k ) x ( k | k ) . 

The next derivation shows how the error covariance 
matrix P k is projected from time t k to time t k + j . 

The optimal estimates at times t k and t k+ j are 
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assumed to be related to the parameter vectors by 

x k + l ■ x k + l - x k + 1 ( A I . 5 ) 

*k = x k ^ x k- 

where the tilde designates estimation errors. Subtracting 
equation (AI.2) from (AI.4) leads to 

*k+l * ®k x k * 1 ’k w k • < A I • 6 ) 

The error covariance at time t^+i is 

p k + 1 5 E { Xk+ i^k+ 1 ^ (A1.7) 

= E{(4» k 3c k -r k w k )(a» k x k -r k w k ) T }. 

where the superscript T denotes the transpose operation. 

Using equations (AI.3), and working through some algebra 

T 

results (where P k = E{x k x k )) in the 
covariance update 

P k + 1 = ®k p k®k T * r kQk r k T • (AI.8) 

This is another step in the Kalman filter algorithm, and in 
the notation of Jazwinski (1974) is 

f 

P ( k+ 1 i k ) = 4>(k-l ,k)P(k'k)® T (k+l ,k) + T ( k ) Q ( k ) I’ T ( k ) . 
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How measurements are combined with Kalman filter esti- 


mates to yield updated estimates will now be illustrated. 
It is given that optimal estimates at times immediately 
before ( - ) and immediately after ( + ) a discrete measure- 
ment are related to parameter estimates by the equations 


Sk* 1 -) = x k * &k< + ) 
x k ( ) = x k ■* x k ( - ) 


(AI .9) 

( A I . 10) 


where the tilde again denotes estimation errors. The recur- 
sive estimator is assumed to have the form 


Xk< + > = K 


) * K k y k • 


(AI, 11 ) 


where K'k and Kk are Kalman gain matrices (which 
are yet to be determined). Equation ( A I . 11 ) simply means 
that the optimal parameter estimate (after a discrete 
measurement) is related to some combination of the optimal 
parameter estimate (before the discrete measurement Is 
taken) and the measurement, itself. The relative contribu- 
tion of Xk(-) and yk are determined by the gain 
coefficients K'k and K^. 

Insertion of equations (Al.l), (AI.9) and (AI.10) into 

equation (AI.ll) produces 


&k* + ) " {Kk'-J+KkMk}xk + Kk'Xk( ) + Kk v k* 


(AI . 12 ) 
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Taking the expectation of equation ( AI .12) and assuming 
E { X k ( - ) } = 0 (the estimator derived herein produces 
unbiased results ) , one finds 

E ( x k (+ ) ) = (K k '-I+K k M k }E(x k ) . ( A I .13) 

For the updated estimates to be unbiased, it is required 
that E { x k } =0 ; this can be achieved by setting 

K k '= I-K k M k , (AI . 14) 

where 1 is the identity matrix. Placing equation (AI-.14) 
into equation (AI.ll) yields equation (AI.15) 

x k ( + ) - x k (-) 4 K k ( y k -M k x k ( - ) ) , (AI.15) 

which is another step in the Kalman filter algorithm. In 
conditional probability notation the equation is 

x ( k + 1 j k+ 1 ) = x ( k+ 1 | k ) 4 K(k+1 ) (y k + 1 -M(k + l )x(k+l | k) ) . 

The error covariance can also be updated in a manner 
similar to equation (AI.15). This starts by using equation 
( A I . 9 ) and a variation of equation (AI.15) 

x k ( + ) = { I ~ K k M k } x k ( - ) + K k y k ( A I . 1 6 ) 

to form X k ( 4 ) = x k ( + ) - x k 

= ( I-K k M k )x k (- ) + K k (M k x k +v k )-x k . 
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And x ( + ) = { I -K k M k }x k ( - ) + K k v k . 


( AI . 17 ) 


The updated error covariance is 

P k ( + ) = E{x k ( + )x k ( h) T } . ( A 1 .18) 

Using equations (AI.3), (AI.17) and (AI.18), and the defi- 

T 

nition P k ( ) = E { x k ( - ) x k ( - ) }. one finds 

P k ( + > = (I-K k M k )P k ( - ) ( I -K k M k ) T + 

K k R k K k T - (AI .19) 

The updated error covariance is related to a combination of 
the prior error covariance and the measurement error 
covariance. In the notation employed by Jazwinskl (1970), 
the update equation is 

P ( k + 1 | k + 1 ) = ( I -K ( k+ 1 ) M( k+- 1 ) ) P ( k+ 1 |k) ( I-K(k+l)M(k+l ) ) T 
f K(k + l)R(k + l)K T (kt-l). 

The last, equation to be derived in the Kalman filter 
algorithm is that determining an optimal gain factor K k . 
This is found by minimizing the error in the updated 
covariance P k (+) (actually the sum of the P k (+) 
diagonal elements) with respect to the undetermined gain 
factor K k . The details of this procedure are given 
in Gelb et a 1 , (1974) and yield the result 

K k = P k (~)M k T {M k Pk(- )M k T+ Rk>" 1 . ( A 1 .20) 
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The optimal Kalman filter gain compares the parameter error 
covariance and measurement error covariance. Notice that if 
the measurement covariance is very large (i.e. poor data), 
the Kalman gain will be small, and the updated parameter 
estimate (see equation (AI.15)) will have a value not much 
different from the predicted estimate x ^ ( - ) . Equation 
(A I. 20) in conditional probability notation is 

K ( k+ 1 ) = P ( k + 1 i k ) M T ( k + 1 ) { M ( k + 1 ) P ( k + 1 j k ) M T ( k + 1 ) + R ( k +• 1 ) } 1 . 

All of the Kalman filter (algorithm) steps have now been 
derived. 
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Appendix IT. 


EFFECTS OF FAULT PARAMETERS AND LOCATION ON CHANGES IN 
PRODUCTS OF INERTIA 

For completeness, a short discourse is presented here 
which illustrates the effects of earth quake strike, dip, 
slip, colatitude, longitude and depth on the changes in 
the products of inertia A C j 3 and AC23 . This section 
closely follows the work of Dahlen (1971, 1973). 

Equations V.3 and V.4 (presented previously ) relate 
changes in the products of inertia, AC-^ and AC23, to 
strike, dip, slip, colatitude, longitude and depth. Dahlen' s 
(1971, 1973) orientation conventions for strike, dip and 

slip are as follows: the strike angle a is measured counter- 

clockwise from North (at the Earth's surface); the dip 6 is 
the angle between the Earth's surface and the fault plane; 
and the slip angle \ is measured counter-clockwise from 
horizontal (in the fault plane). 

Dahlen (19 71) calculates a product of inertia a m p 1 i - 
2 2 1/2 

tude AC = (AC j 3 + AC23 ) for several types 

of faults, with various combinations of c 0 latitude, strike 
and dip. Specific test case results (Dahlen, 1971) are 
presented for a vertical, strike-slip fault (e.g. the San 
Andreas fault) at a depth of 20 kilometers (Figure A I I . 1 ) 
and a shallow dip (6 - 20°) thrust fault (e.g. Chilean 
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subduct ion zone) at a depth of 20 kilometers (Figure All. 2) 

It is evident on Figure A 1 1 . 1 (strike slip fault) that 
the AC amplitude is largest for a colatitude of ninety 
degrees, ignoring the effects of strike orientation; i.e. 
a strike slip fault at the equator has the largest relative 
influence on AC when considering the effects of colatitude. 
The fault strike has maximal effects on AC (for a strike 
slip fault) at integer multiples of ninety degrees (a = 0 
degrees, 90 degrees,...). It appears that a North-South 
trending strike slip fault has the same effect on AC as 
an East -West trending (strike-slip) fault. 

The effects of a shallow angle thrust (dip slip X. = 9 0 ° ) 
fault on AC are shown in Figure All. 2 ( Dah 1 en , 1971). 

AC is at a maximum for a colatitude of 60 degrees and for 
a strike angle of 90 degrees. Thus a shallow angle thrust 
fault at 60° colatitude and trending east-west has 
maximal effects on AC. It appears from the Figures 
(AJJ.l and All. 2 - for the cases studied here) that a 
strike -slip earthquake can produce slightly greater 
changes in AC than the shallow angle thrust fault; one 
may notice that longitude is not varied in Dahlen’s 
(1971) Figures. Knowledge of the symmetry of the earth 
would dictate that AC is not strongly dependent on 
longi tude . 

Up until this point, little has been said about the 
dependence of AC on the depth of an earthquake. AC is 
dependent on depth h through the canonical functions, T j ( h ) 
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Figure All. 2. plot showing variation of AC 
(ACJ 3 * AC 23 ) with 

respect to fault strike a and co latitude 0 O 
for a shallow angle (6 =■ 20 )thrust fault 
of unit slip-area (from Dahlen, 1971) 
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r 2 ( h ) and r 3 ( h ) . The canonical functions are displayed 
versus depth h in Figure All. 3 (Dahlen, 1973). The 
canonical functions all increase with increasing depth, 
with I' 1 and I ' 3 increasing much more linearly than 
1’ 3 . The simplest conclusion that can be drawn is 
(ignoring the other parameters) that AC increases with 
depth h . 

However, it should be remarked that subducted 
(typically dip slip) lithosphere can extend to far 
greater depths (« 700 km) than lithospheric plates 
typically participating in strike-slip motion («100 km). 
On the basis of this observation, it would be valid to 
conclude that AC would be much larger for subduction- 
related earthquakes than for strike -slip events. 
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Figure All. 3. Plot of I’i(h),-r 2 (h),r 3 (h) 

for SNREI Earth model 8073AW. (from I)ahlen.l973 

for details about 8073AW, see Oahleu, 1973) 
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